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I. INTRODUCTION 


A. PURPOSE 


Since the end of the Cold War, there has been a paradigm shift in detecting, 
tracking, classifying and ultimately neutralizing subsurface threats. In its most general 
form, the US Navy now sees a more probable warfare scenario in coastal waters as 
opposed to the Blue Water confrontations occurring earlier last century. The major threat 
of the once large Soviet Nuclear Submarine Force has now tumed into a small-country 
coastal diese] submarine threat. Indeed, the very name for this warfare area changing 
from Anti-submarine Warfare to Undersea Warfare implies a change in our subsurface 
wartare strategy. 

With the shift in subsurface threats come different challenges in detection and 
classification. No longer can acoustic detection systems exploit deep sound channels on 
a relatively noisy nuclear submarine threat. In this new era where the threat of the once 
large nuclear fleet from the Former Soviet Union is limited at best, a new foe has 
emerged in coastal waters. Not only do Naval USW systems have to contend with 
increased challenges such as bottom bounce and added noise due to a higher density of 
neutral shipping, but the platforms of choice are the quiet diese] submarines while on 
battery power. In addition to this diesel threat in a Littoral Warfare environment, mines 
have become more of a threat. Because of their high destruction-to-cost ratio, laying 
mines to thwart US global interests is a cost-effective way to inflict maximum damage 
with little cost. Indeed, the US Navy has suffered more battle damage to their ships due 
to mines than any other subsurface threat since the Korean War. Mines are cheap, easy to 
deploy and hard to detect. To many countries, mines are a viable substitute for an 
otherwise weak coastal defense. 

The nature of Littoral Warfare possesses unique challenges to detecting and 
classifying subsurface threats. Since the Walker Treason in the middle 1980's, enemy 


acoustic signatures have been suppressed to an all-time low. One way to regain a 


considerable edge in Undersea Warfare is to investigate more robust methods of detecting 


enemy acoustic signatures within heavy background noise. 


B. GOALS 


This thesis 1s part of a larger project that addresses the issue of improving the 
tools available to USW systems operators of surface and subsurface platforms. 
Currently, USW systems operators use variations of a spectrogram as their main visual 
display to monitor passive acoustic signatures. This type of analysis is good for signals 
that are long in duration. Existing displays show a long-term history of acoustic data, 
which is good for determining underwater activity over a long period of time. Using the 
same techniques, short-duration signals (or transients), however, may go unnoticed. For 
example, if an operator is either fatigued or focusing on long-term trends, he might 
overlook a very slight indication that a transient signal a fraction of a second in duration 
had just occurred. 

The thesis goal is to create more robust visual and audible features at the 
operator’s disposal to enhance the detection of transient signatures. Specifically, the 
thesis will develop and investigate three signal processing procedures based on Wavelet 
and Multiresolution analyses, which will attempt to detect transient signals currently 
overlooked by common signal processing techniques. 

This thesis will also investigate an automation procedure, which will alarm the 


operator of such transient signals. 


C. METHODOLOGY 


It is the main thrust of this thesis to detect transient signals in a noisy acoustic 
environment. One way to investigate the improvement of detecting desired acoustic 
signatures imbedded within a noisy environment is by studying the statistics and nature 
of the signals. Although littoral warfare might be difficult in terms of background 


shipping and silent threats, one positive aspect is that acoustic data is easily accessible 


to 


because the surveillance area is small and there exist well-known shipping lanes and 
underwater topographies. 

Current methods of transient detection use time-frequency analysis and stochastic 
modeling. Although these methods have been proven quite useful, they have their 
limitations. This thesis will investigate the root causes of these weaknesses and explore 
some ideas to circumvent them. Ideas such as Filter Banks using Wavelet Transforms 
and Multiresolution Filtering just might be the missing tools to enhance transient 
detection. 

Three methods have been developed to detect transient signals in noise. One of 
these methods is based on Wavelet Analysis and the other two based on an ARMA model 
and multiresolution processing. A standard data set is used to compare the new signal 
processing methods against well-established methods. From this data, conclusions will 
be drawn and recommendations will be made for further improvement of acoustic 


transient detection systems. 


D. BENEFITS 


The confidence of a decision is directly proportional to the amount of concurnng 
evidence. Consider the method of detecting submarines using the classic method of 
Target Motion Analysis (TMA). Since TMA uses nothing but passive techniques, 
determining a target’s bearing, range, course and speed can be challenging. What 
strengthens this method of detecting underwater weapons systems is the redundancy of 
evidence. Information that does not show up on a Dead Reckoning Trace (DRT) plot 
might be found on the time-frequency table or on a Moboard plot. Given a well-behaved 
target, the success of TMA relies on the fact that dynamic changes will be detected in 
different domains. For example, a shift in frequency may be the first detection of a 
change in course or speed, but probably won’t be the first method to estimate a target’s 
range or course. The different techniques in TMA are independent in the sense that they 
do not rely on similar assumptions. These methods each provide a piece of a puzzle, 


which when put together, can accurately determine particular warfare scenanio. 


Like the set of tools used in TMA, an acoustic transient detection system needs to 
be developed which takes the same fact-layering approach in order to make an informed 
decision. This thesis does not attempt to improve on well-established signal processing 
techniques such as time-frequency analysis. Rather, the goal of this thesis is to create 
additional layers (based on fundamentally different assumptions) and add to the fact- 
finding set in order to enhance the robustness of detecting transient signals. 

The benefits of this thesis will be the addition of new tools, which will enhance a 
transient detection system. One of the new signal processing techniques may be the 
necessary tool needed to detect some transient signals. 

Another benefit to the undersea warfighting capability of a ship is to discuss 


methods to enhance the automation of detecting valid transient signals. 


He CURRENT TRANSIENT SIGNAL DETECTION TOOLS 


Although current methods of detection can successfully identify transient signals, 
there is a threshold floor or frequency bandwidth to which signals become 
indistinguishable from background noise. Of those cases where current detection 
schemes miss a signal, we will show that the new methods discussed in this thesis will 
have a strong indication. This chapter establishes strengths and weaknesses of existing 
signal detection methods both in the time and frequency domains. The intent is not to 
produce an all-inclusive list of current signal processing techniques, but rather to 
compare new methods against existing ones. In this way, we hope to show that these new 
techniques will not replace current ones, but rather support them as they add to the 


robustness of the whole Acoustic Transient Detection System. 


A. ESTABLISH A BENCHMARK SIGNAL 


In order to meet the goal of enhancing transient signal detection, a comparison 
between established methods and new methods presented in this thesis must be 
performed. To accomplish this evaluation, a single representative signal will be used as a 
benchmark. The model signal is the sound of a fan blowing for six seconds sampled at 
11.025 kHz. Starting after the third second, a series of 11 tapping sounds (each 
approximately .1 seconds in duration) will occur. These taps will be part of the analysis 
on transient resolution. The first three seconds will be used to study the effectiveness of 
all transient signal detection techniques by varying the amplitude, duration, and 
frequency content of an added transient signal. Figure II-1 shows a block diagram of the 
filtering process of the signal to detect the transient noise. Figure II-2 is a pictonal 
representation of the benchmark comparison setup based on the filtering process 


described in Figure II-1. 
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Figure []-2 Whitened and Filtered Acoustic Signal Benchmark (z[n]) 


1. What Kind of Noise? 


It is widely held that since acoustic ambient noise is comprised of many random 
sources, it can be assumed to be Gaussian [2]. This assumption is based on the Central 
Limit Theorem, which states that as the number of random variables in a sequence 
becomes large, the probability distribution of the coi of the sequence approaches a 
Gaussian random variable [4]. In the case of shallow water, however, ambient ocean 
noise coming from the Infrasonic Band and the Low Sonic Band may dominate, and 
therefore may compromise the previous assumption. 

The Infrasonic Band comprises frequencies between 1 Hz and 20 Hz. This is the 
frequency band that contains blade rate fundamental frequencies of power-driven vessels. 
The Low Sonic Band has a frequency range between 20 Hz and 200 Hz, which is mostly 
caused by distant shipping, especially in the frequencies around 100Hz [2]. 

At any rate, noise can be defined simply as any part a signal that you don’t want. 
For the specific case of detecting transient signals in shallow water, it is a good 
assumption that the above-mentioned ambient noise bands dominate the noise model. 

Take for example the scenario of attempting to detect a small ship in shallow 
water laying mines among a fleet of fishing vessels. It is desirable to detect a certain 
transient signal like a splash or machinery among many acoustic sources with similar 
frequency content. Since we know that mines are deployed in heavily traveled 
waterways, you can assume that the background noise (defined as the part of the signal 
you don’t want) will have a dominant low frequency content. The model therefore, will 
not be white noise, but colored noise since it has uneven frequency content. 

The colored noise assumption is reflected in the Benchmark Signal. Much like a 
power-driven vessel, a fan blowing produces a similar low frequency content due to the 


blade rate of the fan as it rotates like a propeller. 


2. Are Shallow Water Acoustic Signals Stationary? 


Like most random processes occurring in nature, the background ocean ambient 
noise 1s not stationary. In the example above, one noise source out of many may change 
its course or speed. So by one slight change in frequency, a change in the stochastic 
nature of the ambient noise will take place. However, as in many applications to linear 
filtering techniques, a process can be assumed stationary over a short period of time [5]. 

The benchmark signal for this thesis assumes that the signal is stationary for the 
entire 6 seconds. As the results will show (Figure II-2), assuming the data set is 
Stationary over this time frame is sufficient to distinguish the transient signal from the 
background noise. Like the Benchmark Signal, the short-time stationarity assumption for 


underwater acoustics will apply. 


B. STOCHASTIC SIGNAL MODEL 


The first of two current methods to be discussed is an analysis in the time domain. 
We wish to transform the input signal through a filter to extract a transient from noise. 
To do this, we rely on the innovation process where the output of the filtered sequence is 
white noise. 

It is well known that a process whose complex spectral density function satisfies 
the Paley-Wiener Condition can be modeled as an Autoregressive-Moving Average 
(ARMA) process [5]. One of the properties of the signal model is that it can be reversed 


such that it becomes a whitening filter. Furthermore, since its transfer function H__,(z) is 


fant : : = : , : 
minimum-phase, its inverse H. (z) will be causal as well. Figure II-3 shows a diagram 


of the innovation representation. 


White noise 


x[n} 


White noise 


(innovation process) 





Figure IJ-3 Innovation Representation of a Random Process. (a) Signal Model 


(b) Inverse Filter [5] 


Ih Autoregressive — Moving Average (ARMA) Model 


Any regular stationary random process can be represented as the output of a linear 
shift-invariant filter driven by white noise [5]. Starting with the assumption that our 
input signal x[n] in Figure II-3(b) is colored noise, we want to design a filter H(z) such 
that the result will be white noise. The innovation process will whiten the tme domain 
signal except where the transients are located because it is not part of the innovation 
process. A second look at Figures II-1 and I-2 will show this process. 

We use an autoregressive moving average (ARMA) model in order to represent a 
random process. The ARMA model involves both a nontrivial numerator and 
denominator and when compared to the AR or MA models alone, they often require 


considerably fewer parameters to match a desired power spectral density function [5]. 


Without a priori knowledge of the autocorrelation function, we estimate it 
directly from the data. The Autoregressive model] estimation can then be wnitten as a 


weighted sum of prior values of the observations. Written as an equation, we get 
Minh ated — a, xin o> ein — P|. () 
where the error is 


e{n] = x[n]-—x[n] = y aed: (2.2) 


In equation 2.1, %[n] is the estimate of the n” element. As in many applications, 


we want to minimize the error by applying the Orthogonality Principle which states that a 


vector of weighting coefficients a can be chosen to minimize the mean-square error 
E{|x — 3°} if the error 1s orthogonal to the observations. The final result is the least 


squares form of the Yule-Walker equations, 


(x Xa = A (2.3) 


where 


x{0] 0 ve 0 
x{1] x{0] vee 0 
sed x{P-l] -- x{0] 
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0 aN Ue ed ae] 
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and Sis the sum of the squared errors. It follows then that an estimate for the 


prediction error variance is 
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where J, is the size of the data set. 


The observation matnx X is set up such that it assumes the observations 
outside of N, data points are zero, which is analogous to performing a rectangular 
window operation on the data set. Although this might not be the best method to use in 
terms of pole placement of the vector a, the matrix structure has a nice property since 
(X"7X)is Toeplitz. This property means that the matrix will be positive definite, which 
guarantees that the whitening filter is minimum phase and therefore is a causal stable 
system with a causal stable inverse. This property is important because we are interested 
in the innovation process. Finally, using the autocorrelation method of determining the 
Toeplitz correlation matrix (X”’ X), fast algorithms for determining the elements of a 
can be used such as the Levinson Recursion [5]. From the AR coefficients found from 
equation 2.4, the Moving Average coefficients b can then be computed by any number of 
algorithms such as the Basic Prony Method [5]. The ARMA routine used in this thesis 
uses an algorithm to determine an optimum number of poles and zeros, P and Q, 
respectively. 

As an example, take the Benchmark Signal and find the AR filter coefficients a 


and the MA coefficients b where 


a, 
Ee anGen: = 


1] 


The causal ARMA filter is then 


2} ae a2 
Dat D. 2 eda tO 


FL arma (z~)= (ae); 
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To realize the time domain solution to the problem described in Figure II-3(a), we 


Q P 
xt Sb, w[n — q]- yy eee il aye | (2.6) 
q=0 p=l 


Inverting the model as in Figure II-3(b) the causal inverse ARMA filter is simply 


I+a,2° +a,27 +--+a,2" 


H. (20) =————_——-+—_.. ea 
aru | by + bz +b,27 +--+ byz ; 
The time domain realization to the innovation process 1s finally, 
P Q 
w[n] = S'a,x[n = \|= ¥' bwin apie (2.8) 
p=l g=0 


Running the Benchmark Signal through the inverse ARMA filter (as depicted in 
the model in Figure II-1) results in the innovation process, which is white except where 
the transients occur in time. This is because the ARMA model seeks only to model the 
colored noise as the innovation process. Since the colored noise input has a stationary 
localized frequency content, the ARMA model can do a good job at placing poles to 
model the signal. All other signal content that is statistically dissimilar to the colored 
noise 1s rejected and appears as a “spike”. 

As acheck, in listening to the fan blow, a low frequency “hum” is distinguishable 


and the Il tapping sounds are also recognizable. After the innovation process, the 


background noise sounds like a heavy rain hitting the ground. 


When plotted the 


transients, indistinguishable in the time domain before the innovation process, are now 


quite visible against the white noise. Figure IJ-4 shows the Benchmark Signal before and 


after the innovation process. 
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Figure I]-4 Innovation Process of Benchmark Signal (a) Input (b) Output 


2: Strengths of the ARMA Model 


Since the innovation process 1s based on a statistical model, the frequency content 
of the signal does not affect its detection, even if the signal contains the same frequencies 
as the noise. The ARMA model provides a representation of the noise and it can detect 
signals that are not compatible with the model. The ARMA method picked up on all the 
transients and the height of the spikes of the output signal is proportional to the 
“loudness” of the audible signal, which might be an attractive feature to compare the 
magnitude of two signals. This technique seems to perform well on very short transients 
and it can resolve signals close to each other. Notice, for example, the cluster of signals 
just before the third second in Figure II-4 (b), where the ARMA model was actually able 
to distinguish between three transients less than one-tenth of a second apart from each 


other. 


3. Weaknesses of the ARMA Model 


The reverse ARMA (innovation) process 1s a procedure that produces white noise 
from a colored noise input. This process, as we shall see, 1s very effective in detecting 
transient signals that are not part of the white noise model, appearing as spikes above the 
noise threshold as seen in Figure II-4 (b). However this method, although performing 
well on both narrowband and broadband signals, has its limitations when the power of the 
transient signal is at or below the modeled noise power. As we shall see, the fundamental 
idea of the methods presented in this thesis is to lower the white noise power of the 


innovation process in order to extract transients with smaller signal power. 


C. TIME-FREQUENCY ANALYSIS 


The second method to be discussed introduces a technique that combines both 
time and frequency analysis. The process is like a sliding window in ume taking a 
snapshot of the frequency content of a signal, assuming that the statistical properties do 
not change within the time duration. To view the results, frequency is plotted as a 
function of time and is called the spectrogram. For USW, a similar technique is used 
called the “lofargram,’ which is currently the standard method for sonar signal 


processing [5]. 


1. Spectrogram 


The Spectrogram is based on the discrete-time Short-Time Founer Transform, 


which is an extension to the basic Discrete-Time Founer Transform and is defined as 


A) ¥ x(n) w(n —m)e!" , (2.9) 


imN==-co 


where 
w(n —m) is the sliding window function and 


X(n,w) = ¥ x(m)e?" is the DTFT. 


m=-co 


Completely analogous to the DTFT, the discrete-time STFT is continuous in 
frequency so for digital processing, a method to discretize frequency content like the FFT 
can be implemented to make this method computationally practical. Figure II-5 shows a 
diagram of the STFT process. Along the horizontal axis is a sliding window process as a 


function of time. At each 7,.. +Ar, a “snapshot” of the frequency content of the signal 


at that instant in time is taken. Another way to visualize the process is to look at the 


frequency domain. Along the vertical axis is a set of filter banks, each of equal 


bandwidth Af. The signal is passed through the filter bank at each 7, + At and the 
spectral content within the particular bandpass filter 1s calculated. The latter pictorial 
representation of the STFT 1s particularly useful in drawing an analogy to the Continuous 


Wavelet Transform discussed next. 
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Figure II-5 Time-Frequency Plane Corresponding to the Short-Time Fourier 


Transform [7] 


The (discrete) STFT on the Benchmark signal is provided in Figure II-6. The 
darker colors on the spectrogram represent frequencies with high power. Notice the low 
frequency content below 300 Hz is dominant throughout the entire 6-second sample. 
Constant frequencies at 600 Hz and 900 Hz are also noticeable. The most revealing 
information on the spectrogram is the vertical lines in the latter half of the signal 
indicating the presence of transients. What’s more, the discrete STFI can provide 
information on the type of transient it is. In this case, the transients are broadband 
signals. In addition to the standard Benchmark Signal, three narrowband transients were 
added to the first second to study the effects of narrowband transient resolution. The 
spectrogram located two of the three signals, one at around 5 kHz +500 Hz and one 
around 2 kHz + 500 kHz. The third narrowband transient, however, was not detected 


because it was masked by the low frequency noise. 
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Figure II-6 Spectrogram of Benchmark Signal 


ae Strengths of the Spectrogram 


Among the many strengths of the spectrogram is its colorful visual representation. 
On a single plot, it provides information in both the time and frequency domains. This 
representation therefore gives not only the time at which the transient signal occurred, but 
also its frequency content, which may be vital in classifying the source. For applications 
in sonar signal processing, the spectrogram is commonly referred to as the “lofargram’”, 
which provides a useful “fingerprint” mapping of unique acoustic sources [5]. From this 


display alone, useful information such as detecting dominant narrowband _tonals, 





harmonics and blade rates can be seen. Seeing the frequency content of a contact change 
over time can also indicate a dynamic event such as a change in course or speed. 

The spectrogram works best on broadband transient signals because there is a 
larger contrast of darker colors in a vertical line than there is in the single point 
representation of a narrowband signal. In addition, we have shown that it is possible to 
have a narrowband transient signal masked by the frequency content of the background 
noise. Unless the background noise is completely white, at least some of the frequency 
content of a broadband noise is sure to appear. 

Compared to the ARMA model, the STFT only assumes stationanty over the 
duration of the sliding window only. Recall that the ARMA process attempts to model 
the signal through the statistics of the entire 6-second data set. If it were to assume a 
shorter data set to compute the autocorrelation matnx, the whitening model will be 
different and may or may not provide better results. A further discussion on the window 


size will be discussed in the final chapter. 


i Weaknesses of the Spectrogram 


Like the innovation process, the discrete STFT must assume stationarity for the 
duration of the windowing function. There exist natural and man-made signals for which 
its spectral content is changing so rapidly that finding a proper windowing function is 
difficult because there may not be any obtainable time interval for which the signal is 
stationary [6]. The result is a transient signal may go completely unnoticed. 

Even if it is a good assumption that a short time duration is stationary, there is a 
tradeoff to the resolution in the frequency domain and vice versa. This is called the 
Uncertainty Principle or Heisenberg’s Inequality, which states that resolution in time and 


frequency is bound by the inequality 
] 9 
Ar Af 2—., (2.10) 
470 


with Ar and Af the duration of the signal in time and frequency respectively. 


19 


If localization in time is necessary, Atcan be made small but according to the 


inequality (2.10)Af will become greater and therefore the spectrogram will lose 


frequency resolution proportionately. 

To further illustrate this principle, the two visible narrowband transient signals in 
Figure II-5 are about one-tenth of a second in duration and contain only one frequency 
each, 5 kHz and 2 kHz respectively. From the figure, both transients are stretched 
vertically, blurring the actual frequency spectrum while the time duration at which they 
occur remains short. If a broad window were used in the discrete STFT, the frequency 
would resolve to a single point while the time of the signal would be stretched and 
distorted. 

The spectrogram merely shows the frequency of a signal during a time duration. 
Since the spectrogram is not based on statistical modeling, it cannot “extract” a signal 
from the noise and therefore is unable to distinguish between signal and noise if they both 
have the same frequency content. This 1s particularly important in cases where an 
acoustic transient signal is masked by the colored noise. Consider the Benchmark Signal 
Spectrogram in Figure II-5. Also included in the signal 1s a sine wave at a frequency of 
100 Hz. Since the power of the colored noise is highly concentrated at frequencies below 
300 Hz, the transient signal goes undetected. Unmasking transient signals in this region 
become important for sonar operations because much underwater activity may take place. 

Finally, while the innovation process forces a signal into a white noise process 1n 
order to detect transients, the spectrogram takes advantage of the signal’s frequency 
content over time. If the input signal is white noise, nothing 1s distinguishable since the 


frequency content is spread evenly over all frequencies for all time. 


Hil, WAVELET ANALYSIS 


Since we know that the STFT is resolution-limited in both time and frequency 
simultaneously, we wish to investigate other means to localize transient signals in time 
while somehow still being bounded by the uncertainty principle. In this section, we start 
with a discussion of the Continuous Wavelet Transform (CWT) as it relates to the Short- 
Time Fourier Transform. From there, we discuss the development of the Discrete 
Wavelet Transform, its usefulness with filter banks, and the idea of multiresolution 


analysis, which is the basis of the methods of detecting transient signals. 


A. CONTINUOUS WAVELET TRANSFORM 


To describe the Continuous Wavelet Transform, we start with the equation for the 


Short-Time Fourier Transform 


X(f)= x(t)w' (t- Te?! dr. (3.1) 


—oo 


The integral can be viewed as an inner product operating on data inside the 


windowing function, which measures the “similarity” between the function x(t) anda 


sinusoidal basis function. The result turns out to be the energy level at each specific 
frequency. STFT analysis is useful for narrowband periodic signals because of its ability 
to detect strong similarity between the input signal and the sinusoidal basis function at a 
particular frequency, but it has its limitations to the larger family of non-stationary, 
broadband signals. Based on Fourier Transform restrictions, it 1s desirable to develop a 
different transform that can overcome the resolution limitation. 

In order to investigate a different set of transforms that exhibit these properties, 


we need to look at other possible basis functions. To begin, assume that the width of the 


frequency resolution of the transform // is no longer constant but rather the ratio of Af to 


the frequency f , resulting in 
A 
Af =c. (3.2) 


This constant ratio means that as the analysis of frequency content of a signal 
increases, So does the width of the frequency resolution. Looking back at the diagram of 
the time-frequency plane of the STFT (Figure II-5), the vertical axis represents a set of 
constant-width bandpass filters in a filter bank. In the case of the CWT, the constant- 
width frequency resolution is replaced by a filter bank of constant relative bandwidth 
(called “constant-Q” analysis). This process, in essence, allows the resolution in At and 


Af to vary in the time-frequency plane in order to obtain a multiresolution analysis, 


which will be discussed in greater detail in the next section. Another way to look at the 
filtering process is that the filter bank used is not linearly spaced along the frequency axis 
(as for the STFT case) but spaced evenly over a logarithmic scale [7]. Figure II.1 shows 


this process. Now substituting (3.2) into the Heisenburg Inequality (2.10) gives 


Ni (3.3) 


By looking at equations (3.2) and (3.3), one can see that as the frequency becomes 
greater, the resolution in time decreases and the resolution in frequency increases. 
Conversely, as the frequency becomes smaller, the resolution in time becomes larger and 
the resolution in frequency becomes smaller. This is a fundamental property of the CWT. 
The time resolution of the CWT becomes arbitranly good at high frequencies while the 


frequency resolution becomes arbitrarily good at low frequencies. 


to 
to 


Constant Bandwidth (STFT case) 
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Figure III-1 Filter Division in the Frequency Domain. (a) Uniform Coverage (STFT) 
(b) Logarithmic Coverage (CWT) [7] 


The definition of the Continuous Wavelet Transform is thus 
CWT, (1,4) = | x(t)h’ ax (dat, (3.4) 


where 


Heli) “an (S (3.5) 


vel 


is the basis function called a wavelet. 
A description of a wavelet function will be addressed in the next section. For 


now, some of the properties of the transform will be discussed. From the definition, it 
can be seen that ais the “scaling” factor described in equation (3.2) and 1/ /|a| is used 


for energy normalization. The numerator of the wavelet function shows the time — shift 
property of the sliding window function. 

By way of analogy, the STFT is plotted as frequency vs. time as a spectrogram. 
The CWT is a function of time and scale and therefore will be plotted as scale vs. time as 
a “‘scalogram”. Figure III-1 is a comparison of a CWT to a STFT of the same signal. 


Notice that the CWT of the delta function converges to a single point at t =¢, and the 


energy 1s spread out with increasing scale. The STFT of the delta function shows the 
resolution limitation in the time domain with equal energy at all frequencies. For the 


sinusoidal signal of three frequencies at f,, 2f, and 4/,, the STFT shows the constant 


frequency in time, bound by the resolution limitation in the frequency domain. For the 
CWT, however, we see the effect of the logarithmic filter bank where the low frequency 
resolution is narrowest and continues to double in size as the frequency increases. Figure 
If{-2 shows clearly the principle that the time resolution of the CWT becomes arbitrarily 


good at high scale while the frequency resolution becomes arbitrarily good at low scale. 


a) CWT of Delta Function b) STFT of Delta Function 
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c) CWT of 3 Sinusoids 
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Figure III-2 Comparison of the CWT to the STFT [7] 


B. DISCRETE WAVELET TRANSFORM 


At this point a discussion of the properties of the Continuous Wavelet Transform 
has been made as it relates to the STFT and the benefits in localization in both time and 
scale have been presented. The CWT also reduces the number of coefficients necessary 
to define the signal. During the discussion, ideas of orthogonality, energy distnbution, 
logarithmic filter banks, scaling functions, and wavelet basis functions emerged as 
important properties of the CWT. This section will continue to build on these concepts as 
they apply to the Discrete Wavelet Transform (DWT). 

For reasons of clear analysis and better understanding, we begin with the 
definition of a linear decomposition. A real-valued function can be expressed as a linear 


decomposition such that 


f= Maw. (3.6) 


Parameter a, is the real-valued expansion coefficients, and y(t) is a set of real- 


valued functions of ¢ called the expansion set. When the expansion is unique, the set of 
functions y,(t) is called a basis and it is used to represent f(t). If the basis is 


orthogonal, 1.e. 
(W..%)= [vO Mat =0, k#l (3.7) 
then the coefficients can be determined by the inner product 
a, =(f().W,.) =| fOv, Mar. (3.8) 


For the wavelet expansion, however, a two-parameter system is used such that the 


basis used to represent f(t) 1s now 
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FO= Da ia (1), (3.9) 


where 


j represents the integer scaling index of the scaling function, 
kK represents the integer time translation, 


y(t) are the wavelet expansion functions that usually form an, 


orthogonal basis, and 


a,, are the set of expansion coefficients. 


The wavelet basis set y,,(¢) can be defined in terms of scaled and translated 


versions of the Mother Wavelet y(r), i.e. 


i 
y(t) = 22y(2/t-k). (3.10) 


tw fm. 


Here, the term 2° is a normalizing factor as j increases. This property described 
in equation (3.10) is called the multiresolution condition, which states that if a set of 


signals can be represented by a weighted sum of g(t —k), then a larger set (including the 
original) can be represented by a weighted sum of @(2t—k). In other words, “if the 


basic expansion signals are made half as wide and translated in steps half as wide, they 
will represent a larger class of signals exactly or give a better approximation of any 


signal” [9]. 


i Signal Spaces and Multiresolution Analysis 


To discuss the DWT, it is best to start with the idea of resolution to define the 


effects of changing scale. To that end, we start with a scaling function g(t), which then 
will define the wavelet basis function w(t). We define a set of scaling functions in terms 


of integer translates of the basic scaling function by 


Q,(t)=(t—-k) ke Zand geL. (San) 


where the L” space is defined as the set of all functions f(t) with finite energy, 


+00 


[[f@/' ar <e. 
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The subspace spanned by the functions @, (t) 1s V, such that 
v, = Spanig, (t), ke Z}. (3.12) 


That is to say, the subspace v, 1s spanned by the closed set of all functions 9g, (¢) 
over all integer values k. From equation (3.12) we can say that a function f(t) within 
the subspace v, can be expressed as a linear combination of the set of functions @, (ft) 


such that 


f= a.9,(.)  forany f(rje vy. (3.13) 


From here, we consider a two-dimensional family of functions and derive a similar 


linear combination as in equation (3.13). We start with the definition of a scaling 


function parameterized by scaled and translated versions of Q(t) similar to equation 


(3.10), resulting in 


Z 
2 


Gi C=2792 tk), kEeZ (3.14) 
where 
v, = Spanig, ,(t), k€ Zf. (3.15) 


In the same way Vv, is expressed as a linear combination of the set of functions 


Y,(t) in equation (3.13), f(¢)€ v, can be expressed as 
ft) = Ma,p(2/t +k). (3.16) 
k 


From here we turn to a useful property of the wavelet expansion set [9]. There is 
a basic constraint of multiresolution analysis, which requires a “nesting” of the spanned 


spaces as 
Ge men ei uie yan.) 
or simply 


I 


forall je Z. (3.17) 


By this observation, in order to ensure that the elements in v, are scaled versions 


of the next higher space v,,, we add the constraint, 
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Ge vs = 2) © Ve (3.18) 


This is done to ensure that analysis and synthesis can be performed. That is to 
say, if we conduct a wavelet decomposition at a certain resolution (or scale), there exists 
a reversing process whereby a signal can be reconstructed (synthesized) to its original 
state. Figure III-3 shows a Venn diagram of the nested vector space property of 


multiresolution analysis. 


I? D+: Dv,Dv,Dv,D Vv, 
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Figure III-3 Nested Vector Spaces Spanned by the Scaling Functions [8] 


This set of nested vector spaces now meets the multiresolution condition, which 


states that if a set of signals can be represented by a weighted sum of @,,(7), then a 
larger set can be represented by a weighted sum of g,,,,(¢). In order to relax notation, 


we dropk and assume that g (7) is dependent on the time translation. The nesting of the 
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spans as shown in Figure III-3 and by equations (3.17) and (3.18) is obtained by requiring 


that @(t) is not only in the space v, but also in v,, which is spanned by the function 
~, (t). From equation (3.16), we can see that that the function @,(t) within the subspace 


V, can be expressed as a linear combination of the set of functions @, (t) such that 
g(t) = ¥a,9,(t). (3.19) 


Substituting @,(¢) and @, (t) using equation (3.14) we obtain an equation in terms 


of scaled and translated versions of the scaling function as 


g(t) = ¥ h(n)V29(2t-n), ne Z. (3.20) 


This equation shows a telescoping structure of the scaling function and is called 


the Multiresolution Analysis (MRA) equation, where h(n) represents the weighted 


coefficients called the scaling function coefficients. 
We now have a recursive scaling function spanning a subspace, which is also a 
subset of the next higher scale. There is one more feature of a signal that does not use a 


scaling function, but rather uses a different set of functions. called w(t). This set of 


functions is defined as a basis set that spans the differences between the spaces spanned 
by the various scales of the scaling function [9]. As we shall see, these are the set of 
wavelet functions. 

Although wavelets and their corresponding scaling functions are not required to 
be orthogonal, for the purposes of this thesis, we will limit our discussion to the set of 
orthogonal basis functions. The advantages of having an orthogonal basis are that the 
wavelet expansion coefficients are easy to calculate and are governed by Parseval’s 
theorem which allows the partitioning of signal energy in the wavelet transform domain, 


like that of the classic Fourier transform domain. In vector function space notation, we 


Sy) 


define w, as the orthogonal complement of v, in v,,,. That is to say, \w si is the set of 


all functions that span v,,, and are not included in v,. v, is also be considered 


uncorrelated with w, and the union of the two vector spaces span Vv ,,,. Mathematically, 
all functions spanned by v, are orthogonal to w,. Therefore the basis functions g,, and 


y,, are orthogonal as 


(954). j.0)) =] 9;. (Wj, (dt = 0. (3.21) 


Starting at an arbitrary scale spanned by the scaling functions in v,, say 7 =0, 


we Can write 
2 
Vey eV, @ Gala 
Now from the orthogonal relationship between v, and w,, we get 
Ve—=V, oe (3-22) 


where © indicates a union operator. 


Finally, due to the “telescoping” nature of the subspaces v,, j © Z and the 


orthogonality of the corresponding w,, jé Z, all of L* can be expressed as 
L’ =v, @w, Ow, Ow, ©---. ews) 


Figure III-4 shows the Venn diagram of the entire vector space spanned by L’. 
Assuming that our arbitrary starting point in the infinitely telescoping set of nested vector 


spaces spanned by the scaling functions is v,, we can see that w, is orthogonal to Vv, 


and a subset of v ,,, 


1 


as shown in Figure III-5. These wavelets reside in the space spanned 
by the next narrower scaling function and therefore the multiresolution principle applies. 
Since w, CV,, the set of wavelets in w, can be represented as a weighted sum of the 


next higher scale g(2r) resulting in 


with= Yh, (n)V2@(2t —n), We (3.24) 


where h, (7) represents the weighted coefficients that define the wavelet. 


It can be seen that in equation (3.24) the wavelet function is expressed in terms of 
a linear combination of scaled and translated versions of the scaling function, which 
implies a relationship between the scaling function and the wavelet function. Indeed, by 


taking equations (3.20) and (3.24), it can be shown that A(m) and h,(n) are related by [9] 


LE Sie). (3.29) 
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Figure III-4 Scaling Function and Wavelet Vector Spaces [8] 





Figure [fI-5 Orthogonal Basis Set in Subspace v ,,, 


jae Implementation of the Discrete Wavelet Transform 
We now have defined a wavelet system in terms of its scaling function and 
Mother Wavelet. As it turns out, it 1s not necessary to actually know @(t)and w(t). All 


we need are the filter coefficients to the MRA equation in order to perform a discrete 
wavelet decomposition. This section introduces the Discrete Wavelet Transform (DWT) 
and how it can be implemented. 

In the last section we described the whole vector function space as a union of the 


coarsest scale with all wavelet subspaces such that 
L =v, Ow, Ow, Ow, O---. 


Based on this fact, we can now express a function g(t) as a series expansion in 
terms of the scaling function @, (1) and the wavelet functions y,, (¢) giving the general 


form of the equation 


a(t) = Velo. +d Yd, Oy, (0). (3.26) 
k 


k j=j0 


Here, jO is the arbitrary position for the coarsest scale whose space is spanned by 
the scaling function @ ,.(¢).. Notice that the inner summation on the second term begins 
at j0, indicating that the wavelet coefficients begin at the same level as the current scale 


level. If the conditions of orthogonality are met, these wavelet coefficients can 


completely describe the function g(t) and are useful for analysis, denoising, compression 
and perfect reconstruction of the original signal [10]. Assuming that the wavelet system 


is orthogonal, the coefficients c ,(k) and d,(k) can be found by the inner product of @,, 
and y,, respectively. 


Performing the inner product to find c (A) we have 


c,(k) = (20). )40)) =] 309), (dt 


=| g(0)2? g(2/1-k)dr. (3.27) 


[~~ 


Restating the MRA equation, 


eo=> h(n)V29(2t -n), 
we substitute in the next higher scale and translate the function by k, resulting in 
g(2/t-k)= Dho2922"1 —2k—-n)). 
By making a variable change m= 2k +n, we get 
o(2/t-k)= pa h(m — 2k)V2@(2/'1-k). (3.28) 
Now substituting (3.28) into (3.27) it can be shown that 
j+l 


c)(k) = SV h(m—2k)| g(e)2 * G21 -k)dt. (3.29) 


Finally, the integral expression is simply the inner product 
(2(1).Q jem (t)) = Cj), 


and therefore 
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c,(k) = \ h(m—2k)c ,,,(m). (3.30) 
A similar derivation for finding the coefficients d ,(k) yields 
d ,(k) =) h,(m— 2k)c ,,,(m). (3.31) 


Equations (3.30) and (3.31) show that the scaling and wavelet coefficients at the 


next lower scale j can be computed by convolving the scaling coefficients at the current 
scale j+1 with the time-reversed weighting coefficients A(—n) and h,(—n), then 


downsampling the result. Figure HI-5 shows a diagram of this process. 





Figure III-6 Two-Stage, Two-Band Analysis Tree [8] 


3. Analysis in the Frequency Domain 


In the following section, it will be shown that the coefficients of the Haar Wavelet 


filter coefficients h,(n)are a lowpass filtering operation. Likewise the coefficients h, (72) 


are a highpass operation [9]. Together, they form a Quadrature Mirror Filter Bank 
(QMF) system. In addition, a reverse filtering operation can be performed allowing for 
perfect reconstruction of the onginal signal, also known as a synthesis procedure. This 
makes sense because for any transform, a reverse operation must be obtainable for the 
transform to be widely used. 

In the frequency domain, the signal coming in is processed through a filter bank 
that splits the frequency content into two parts. The high frequency band of the signal is 
represented by the wavelet coefficients also known as the details. The low frequency 
portion of the signal is represented by the scaling coefficients also Known as the 
approximations. The lowpass portion can then be processed through the same filter bank 
again. Theoretically, lowpass filtering the scaling coefficients can be repeated 
indefinitely, but in practice is limited to the number of data points in the signal. This 
limitation exists because every time the signal is passed through the filter bank it gets 
decimated by a factor of 2, which effectively means the number of data points in the 
signal is cut in half. 

The first stage splits the frequency band in half, the second stage splits the low 
frequency portion into quarters, and so on. The resulting “analysis tree” makes up a 
logarithmic filter bank as seen in the CWT. A diagram of this process can be seen in 


Figure [JI-7 as it compares to a constant width filter bank. 
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Constant Bandwidth (STFT case) 
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Constant Relative Bandwidth (WT case) 





Figure III-7 (a) Constant Width Filter Bank (b) Frequency Bands for the Analysis 
Tree [8] 


C. HAAR WAVELET 


Up to this point we have studied the Wavelet Transform as a linear 
decomposition. We now move from the abstract to the application. Equations 3.30, 3.31 
and Figure IIJ-6 all describe the process of decomposing a signal into successive details 
and approximations. Now looking at the decomposition as an application to signal 
processing, we see that it is a series of simple filtering and downsampling operations. 
For our particular application, the actual numerical values of the decomposition 
coefficients are not important, only their relative values and so the coefficients can be 
normalized. By viewing the DWT as a simple signal processing procedure, we have 
essentially found a way to split and downsample a signal into mutually independent 
signals. As we shall see, this procedure provides useful properties in detecting transient 


signals by removing background noise. 


ie Haar Wavelet Parameters 


We are now ready to discuss the actual wavelet filter coefficients h(n) and 


h,(n). There are many wavelet basis sets to choose from and a discussion on choosing 


specific types will be discussed in the final chapter. For the purposes of this thesis, we 
will choose the Haar Wavelet system and derive all transient detection methods and 
conclusions based on this set. Admittedly, this is the simplest of all wavelet basis sets but 
the results will show a marked improvement over current detection schemes. In addition, 
it is More intuitive to describe the filtering processes in the upcoming chapters using Haar 
wavelet coefficients. Once the transient detection methods have been discussed, a more 
general case can be easily substituted. 


The Haar Wavelet coefficients are defined as 


h(n) = {h(0),h()} = (3.32) 


Sl- 


l 
AB 
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for the lowpass filter and 


(3.33) 


h, (1) = {h, (0), h, (1) }= I+. 5 


for the highpass filter. 
A more detailed discussion explaining how the Haar Wavelet meets the necessary 


conditions of a Wavelet system is presented in Appendix A. To find the shape of the 


scaling function, recall the MRA equation (3.20) 


g(t) = ¥ h(n)V29(2t - 7). 


Substituting the scaling function coefficients, the scaling function becomes 


g(t) = 9(21) + G(2t -—7n). (3.34) 


Equation (3.34) clearly shows that the Haar Scaling function is a unit-width, unit- 


height pulse function and is equal to the sum of two scaled and shifted versions of itself. 


To find the shape of the Haar Wavelet, we use the wavelet equation (3.24) 


w(t) = » h, (n)V20(21 =f). 


Substituting the Haar Wavelet function coefficients, the above equation becomes 


W(t) = p(2t) — p(2t— 1). (3.35) 


In the same manner as the scaling function, the shape of the Haar Wavelet is the 


sum of two scaled and translated versions of g(t), except now the coefficient on the 
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second term is negative reflecting the negative coefficient in h,(m). A diagram of the 


Haar Scaling and Wavelet functions is shown in Figure III-8. 


Q(t)=9(2t)+ g(2t-1) 


W(t) = 9(2t)-@g(2t-1) 


Figure IT{I-8 (a) Haar Scaling Function (b) Haar Wavelet Function 





For the analysis of both transient detection procedures, it is vital to understand the 
effect of the wavelet decomposition filtering process in the frequency domain. We start 


out by viewing the Haar Scaling filter coefficients in the z - domain as 


Me) Se— Ee) (3.36) 


yy 


She 


Substituting e’” for z we get the frequency response of the FIR filter 


H,(@) = (i+e*). (3.37) 


a 


Likewise, the frequency response of the Haar Wavelet filter coefficients is 


BU) ae (3.38) 


V2 


It is important to note here that the frequency response of the scaling coefficients 
h(n) (also known as h,(n)) is that of a lowpass filter with its cutoff frequency at the 
Nyquist frequency rate (H(z) =0). It can also be shown that the frequency response of 
the wavelet coefficients ,(7) is that of a highpass filter with its maximum at 
H(7)= V2 and its zero magnitude at H(0). This is the property that allows the input 


signal to be filtered into a “constant-Q”’ filter bank resulting in equal highpass and 
lowpass portions preceding every downsample as described in Figure IJJ-7. Figure II-9 


shows the frequency response of the PR filter bank used for the Haar Wavelet system. 





Figure I1I-9 Frequency response of Haar Scaling and Wavelet Functions 


2 Haar Wavelet Summary 


For the DWT, the Haar Wavelet and Scaling functions themselves are only 
important in terms of describing the wavelet decomposition process heuristically. The 
only parameters needed to apply wavelet analysis are the scaling coefficients to the 
Multiresolution Analysis equation. Table [JJ-1 shows a summary of the Haar Wavelet 
system. 

Indeed, the scaling coefficients (mn) and the wavelet coefficients h,(71) work 
together to form a perfect reconstruction filter bank. When a signal is processed through 
the FIR filter /,(7), 1t is being passed through a highpass filter. When the signal 1s 
downsampled, the set of values is the wavelet decomposition coefficients called the 
“details”. Likewise when the same signal passes through the FIR filter h, (7), it is being 


processed through a lowpass filter. After the signal 1s downsampled, the set of values is 


the scaling decomposition coefficients called “approximations”. The approximations can 
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then be decomposed again using the same filter bank and downsampling indefinitely 
(theoretically) to allow for infinite (theoretically) resolution in both time and scale 
(Figure III-6). 

Finally, all details at different scale and the approximation at the smallest scale 
are orthogonal to each other. This property, as we shall see, implies that all signals at the 
output of all filter bank-downsample phases are uncorrelated to each other and is the 


basis for the first transient detection method. 


NAME FUNCTION SHAPE COEFFICIENTS FREQUENCY 
RESPONSE 


] vee 


HAAR (j= 
ot) ‘ otherwise 


SCALING 


1 ee 
2 


HAAR 
WAVELET 


] 
th=<-] -—-St<]l 
y(t) A 


QO otherwise 


} 





Table IJ-1 Haar Wavelet System Summary [15] 
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IV. METHOD 1: INNOVATION - DWT PROCESS 


With the desirable properties of the DWT, we focus on the first method to 
enhance transient detection. This section will present a signal detection design based on 
decomposing the innovation process through the Haar Wavelet system. The desired 
result will be to remove noise through the filtering process while keeping the signal 


intact. 


A. METHOD 1 DERIVATION 


As mentioned in the previous section, the first method takes advantage of the 
orthogonality of subspaces in wavelet decomposition. To say that the details and 
approximations at each stage are orthogonal implies that they are uncorrelated to one 
another. As a proof, consider the whitened input signal processed through one stage of 
the Haar Wavelet system as in Figure IV-1. Note that e,,[/] and e,,.[/] are using the 
same index / because they are produced from the same input and are decimated by the 


same value. 
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Figure IV-1 Wavelet Decomposition on Whitened Signal 


From the previous figure, we can express the approximations e,,[/] and the 


details e,,,[/] in terms of the input signal e[m] such that 


e, p(/] = —(e[2/] + e[2/ -1]) and (4.1) 


No] e 


ell] = —(e[2/] — e{21 -1)). (4.2) 


tWle 


To determine the correlation between the two signals, we perform the cross 


correlation operation, 1.e. 


Efe oll} -€ gol} =~ E((E(20 + ef2 = 1) -tef 21) -€(21-1)} (4.3) 
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Carrying through the multiplication, the expression yields 
| | 2 
E\elidee 0a — 7E| le(20]| JE} le(2/ - 1] iz 0 (4.4) 


Since the signal e[m] is wide sense stationary, E\le(21}° | is equal to 


E ‘le(2I _ vii and equation (4.4) is zero and hence, e,,[/] and e,,[/] are uncorrelated. 


The idea of Method 1 is to perform a DWT on the whitened signal. This setup 
passes a white noise process through the wavelet decomposition where we expect white 
noise at each resolution stage. Based on the properties of the Discrete Wavelet 
Transform and the proof at the beginning of the chapter, all the details and the 
approximations at the same multiresolution stage are uncorrelated. The description here 
is a denoising process where we attempt to lower the noise power by successive filtering 


and downsampling procedures. 


B. METHOD 1 SUMMARY 


Method 1 uses the Haar Scaling and Wavelet filter coefficients for the 
multiresolution filter bank system and is shown in Figure IV-2. Keep in mind this 
method assumes that the DWT is operating on a white noise signal, which means that the 
use of the ARMA model is used pnior to the first stage of the filter bank. If the 
innovation process detects a transient signal, note that there is no need for this algorithm. 
What this method attempts to do is go beyond the detection capabilities of the innovation 
process by removing excess noise in order to enhance the robustness of the detection 


process. 
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Figure IV-2 Method 1: Innovation — Discrete Wavelet Transform Process 


V. METHOD 2: DWT - INNOVATION PROCESS 


We now take a different approach to the same problem using a derivative of the 
Haar Wavelet and its multiresolution capability. As shown in Figure IV-1, Method 1 
took the wavelet decomposition after the innovation process. Now we will investigate a 
method where we perform the multiresolution filter bank process before the innovation 
process as in Figure V-1. The idea of this method is to investigate the effectiveness of 
another denoising scheme. 

Method 2 takes a signal and processes it through the wavelet decomposition. By 
the properties of the DWT, the signals at each channel will be mutually uncorrelated. At 
this point, the decomposed signal is still not whitened. In order to whiten the input, there 
must be an independent ARMA filter for each channel because the same set of ARMA 
filter coefficients is only useful for one channel. Fortunately, it turns out that just as there 
is a relationship between the filters in the wavelet filter bank, there is a corresponding 
relationship between ARMA filters in the innovation filter bank and is dependent on the 
Haar Wavelet filter coefficients. 

This chapter will discuss a procedure to determine a set of whitening filter 
coefficients from the original set of ARMA filter coefficients. The first section starts out 
with modified Haar Wavelet filter bank and begins the procedure to determine all the 


ARMA filter coefficients. 
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A. METHOD 2 DERIVATION 


Figure V-1 shows a block diagram of the multiresolution innovation process. 
First, note that the multiresolution filter bank is a generic highpass/lowpass filter pair that 


are orthogonal to one another. Instead of the Haar Wavelet, the filter coefficients are 


multiplied by a scalar value of _ so that 


/2 


hy (n) = 5 Post (n) = E ;| (Say) 
and 
h,(n) = 5 tn (n) = F 3 Oo) 


The new set of highpass and lowpass filters can be viewed as an averaging 
(lowpass) filter and a difference averaging (highpass) filter. The innovation at the output 
of the whitening filter bank are orthogonal to each other where the “tildes” represent the 
highpass filter channels, or in terms of the DWT, the details. The subscripts represent the 


relative sampling frequency in the sense that e, is decimated by a factor of 2“. The 


dashed lines in Figure V-1 represent a set of cascaded lowpass filters that are 
intermediate steps for further decomposition, yet can also provide useful information. 
Although the decomposition can go on indefinitely, there comes a point where a short 
transient signal is averaged excessively with the background noise. By experimentation, 


three levels of resolution were found to be the most useful in this thesis. 
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Figure V-1 Method 2: Discrete Wavelet Transform — Innovation Process 
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It is also important to note here that the wavelet and ARMA filter coefficients are 
interrelated, and the model is designed exclusively for the Haar wavelet and its 
derivative. A more general form of Method 2 1s possible, allowing the use of any wavelet 
system, but is a subject for further research. 

Consider the ARMA process without multiresolution as in chapter 2. For notation 
simplification, we define the ARMA filter as a transfer function operator on a white noise 


process such that, 


B(z"') 





he e[n], aye 
yol?] es [71] (5.3) 
where 
A(zr ile a, zac coe (5.4) 
JE(C2 stop etd = Doz (5) 


and e[n] is white noise with the second moment being 
E\len]| f=o,7. (5.6) 


From this definition we define y,[/] as the output of the lowpass/downsample 
procedure, also known as the approximation of the DWT. The goal of this section is to 
determine the ARMA filter coefficients of y,[/] as they relate to the ARMA coefficients 
of the orginal signal. Figure V-2 shows a simplified diagram of Method 2 where we 
consider a single stage of the approximation channel. Once the procedure to determine 
the ARMA filter coefficients at this channel is complete, we will go back and consider 


the ARMA model for the detail channel. 
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Looking again at Figure V-2, we want to expand the innovation process in the 
sense that we want to decompose a signal y,[m] as two independent signals as in Figure 
V-2. The transfer function F (27) is the lowpass averaging filter described in equation 


(5.1) so that 


] ] 
y 12] gl Cease (557) 


The sequences eé,[n] and e,[/] in figure V-2 are white noise processes where 
é,[/] is at half the sampling rate as e,[n]. Note here that we define different white noise 


processes throughout the derivation and the terminology can be confusing. The white 


noise signal e,[/] is the final state of the multiresolution filtering process. At this point 
there is no guarantee that e,[n] and e,[/] are uncorrelated and is subject for future 
research. However, it will be shown that the method does provide expected 


improvements in denoising. 





Figure V-2 Single Stage of Multiresolution Innovation Process 
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Starting with the ARMA model of the original signal, we split the original signal 
into even and odd parts. We will define the equation now followed by the proof. The 
goal for this setup is to produce two white noise processes that are independent of one 
another. From the ARMA model in equation (5.3), we wish to define the same process 


decimated into even and odd terms giving 


Cie ane) 


yol2] ]_| Di) DG) | fel at 
Wiel ae Cre) Cle) | le ln 
D(z ) D(z”) 
where 
e,{2] = e[2/], (5.9) 
e [lj = e[2/ —-]]J, (5.10) 
D(z) = A(z" )A(-z"), and ‘exlal 
C(z )\t mG )=B er Az"). Gz) 
where 


C, are the even coefficients of the expression B(z™')A(—z~') and 


C, are the odd coefficients of the expression B(z™')A(—z"'). 


By breaking up the onginal signal y,[2] into even and odd equations, we get an 


identically equivalent system of two equations dnven by two independent white noise 
processes at one-half the sampling rate of the onginal signal. If in effect we can 


construct the sequence y,[n] from two independent white noise processes, then by the 
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properties of the ARMA model described in Chapter 2, the inverse processes will take the 
sequence y,[n] and create two independent white noise sequences. 
The proof for the even/odd decomposition of a signal shown in equation (5.8) 


starts with the ARMA filter determined from the original signal y,[7] such that 


Be) 


H(z"')= me 
A(z ) 


(5.13) 


In order to apply the Noble Identity later on in the proof, all transfer functions will 


=| 


be described in terms of z™ instead of z. We want to represent the ratio of polynomials 


in (5.13) by separating it into even and odd powers of z™' giving us 
Hc) = ae oc (ae (5.14) 


To determine what H,(z"') and H,(z') are, we assume that both transfer 


functions have the same denominator polynomial as a function of z°. We accomplish 


this by multiplying the numerator and denominator by A(—z"'), i.e. 


H(z") ae ) si, eee 
4G AGaD 


Recall that a polynomial can be decomposed into even and odd parts as in 


equation (5.14). Performing this decomposition on the numerator C(z™') we get the 


polyphase decomposition of H(z™')as 


-2 


: — : wees cr) Ch) 
H(z" )=H,(z°)+2°4H, (2°) = + 2 (5.16) 
D(z") D(z~) 








The numerator of equation (5.16) can be thought of as separating all the even and 
odd powers of z™’ from the multiplication of polynomials B(z7') and A(-z7'). To 
show the denominator is a polynomial of powers of z-, however, takes some 
explanation. The denominator is the multiplication between polynomials A(z) and 


A(—z"'). We choose to factor the polynomial A(z7') in terms of its poles Pp; as 


P 


A(z) =|] Q-p2,27), alg) 


i=] 


where the total number of poles is also the filter length P. 


In the same way, A(—z~') can be factored as 
P 
ACz*)=] [G+ pz"). (5.18) 
i=] 


Substituting equations (5.17) and (5.18) into the denominator of the polyphase 


decomposition equation (5.16) we get 


P 


AC iace) — | eeaces) cs 2a) (5.19) 


i=] 


Again, the single series product reflects the fact that the same index is used for 
both expressions. You can see that the cross terms of the expression inside the senes 


product are zero so that the resulting expression yields 


P 
A(z") A(-z") = []- poz") = D(z”). (5.20) 
i=] 
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From the ARMA decomposition equation (5.8) we see that a signal can be 


synthesized from white noise processed through even and odd parts of the ARMA model 
as in Figure V-3. 





Figure V-3 Block Diagram of ARMA Decomposition Equation 


Now consider each output in Figure V-3 separately. Starting with the even terms 
y)(2/] of the original signal we have the block diagram as shown in Figure V-4. Note 
that the downsampling operation is taken to each branch on the other side of the 
summation. At this point we use the Noble Identity as in Figure V-5 and we get the even 


and odd parts of the ARMA model in terms of z™' shown in Figure V-6 [16]. 


7 





Figure V-6 Block Diagram for Even Terms After Applying the Noble Identity 
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From Figure V-6 we see that the white noise process e[n] is split up into its even 


and odd terms as shown in figure V-7. 





Figure V-7 Block Diagram to Determine y,[2/] From a White Noise Process 


The equation to find the even terms of signal y,[m] 1s then 
y,(2]= H(z” )e(21]+ A, (27 )e[2l - 1). (5.21) 
Next we complete the proof by considering the odd terms of the ARMA 


decomposition equation (5.8) as in Figure V-8. Following the same concept used to 


compute the even terms of y,[n], we apply the Noble Identity again leading to the result 


shown in Figure V-9. 
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Figure V-10 Block Diagram to Determine y,(2/ —1] from a White Noise Process 


Therefore, Figure V-10 shows that the expression for the odd terms of signal 


y ln] is given by 
y {22 - = A. (27 )e[21 -1] + 2" H(z" )e[2I]. (5.22) 


Now consider the even and odd terms of the white noise process. Point by point, 
there is no correlation between the two processes and therefore, it can be shown that a set 
of odd and even terms of a single white noise process is completely independent of each 
other. To reflect this uncorrelatedness, we re-define the odd and even parts of a single 


white noise process as two completely independent white noise processes as 


SNPS | AU. (5,23) 


2 a as 9p (5.24) 
It follows then that the complete model of the signal y,[m] is the combination of 


its odd and even parts as defined in equations (5.21) and (5.22) and Figures V-7 and V- 


10. Combining these results along with the newly defined independent white noise 


processes e,[{/] and e,[/], we get the system of equations 


y {(2=H.(z)e (J+ H, (ze, (1) (625) 


y (2i-l =A, (ze, [I] +z A(z” ef, (5.26) 


and a block diagram described in Figure V-11. Finally, by combining equations 


(5.25) and (5.26) into matnx form we get 
vol] )_| H(z") A(z") eal (5.27) 
ee nec ei <> ) | | e,{0) 
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which is the polyphase ARMA of equation (5.8). 





Figure V-11 Block Diagram to Model a Signal from Two Independent White Noise 


Processes 
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We have shown that it is possible to model the even terms of a signal separately 
from the odd using two independent white noise processes. Based on this fact, along 


with the ARMA model given from the original signal y,[n], we wish to determine the 
corresponding ARMA model for y,[/]. Recall y,[/] is the lowpass approximation of the 
original signal y,[n] as described in the multiresolution operation of Figure V-2. 


As it has been said earlier, the multiresolution ARMA filter bank is dependent on 


= ] Z 
the filter F(z~'). So from the averaging filter F(z~ ) = sd +z), y,[/] becomes 


21 
(se = yol21] += yol21- 1) = E Al vol (5.28) 


valet = 


Substituting equation (5.27) into (5.28) we get 


11D He) We) Tet 2 
Fe oe € : ; ao 
si E er “aa ba —_ 


Carrying out the multiplication and substituting H Ree yand H Rae in terms of 


Ce) and Hohe ) 
D(z7) D(z 2 





described in equation (5.16), we get 








eee Cz Coal Vagus) Chee) 
1)=—|— cele 
y (0 Daas a a (I}+ 5 Se SO ben 


Making one final substitution the final expression becomes 








F,(z7') 1 Coe) 
jee 1]+-2 | 
y, [4] Diz “bt Dizry ll (5.30) 
where 
F(z) =—[e,(2)+27C,(27)] and (5.31) 


Ee) =lee" ee Cam )]. (S32) 


Now consider equation (5.30) as a difference equation. We want to determine 


y,[/] as it satisfies the equation 


y= > elit) f, el > dy lm), (5.33) 
n=0 n=0 m=) 


where 


d, = coefficients to the polynomial expression D(z™) assuming d, =1, 


>> 
T 


coefficients to the polynomial expression F,(z™'), 


= coefficients to the polynomial expression F(z‘), 


=> 


= length of original Autoregressive (AR) model A(z”), 

= length of original Moving Average (MA) model B(z'), 
=—i2P =|], and 

BoaO |; 


ES as 
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In order to solve this difference equation, we convert equation (5.30) from a sum 


of two transfer functions into a sum of two state space models. Recall the Type | 


controllable canonical state space form for the Direct Form realization of the transfer 


function [13] 


bz” +b,z" +++-b,, 
iC a,” a ar 
zz t+az te +ay 


palette 0 1 age. v,[n] 
valzt+i]} | 0 Or. “a 0 v, [7] 
ro ee 
Vata — Cn a ee — ly |e 


Therefore, equation (5.35) is of the form 


Vin+l] = A-V[n]+ Bx{n], 


and 


yr] = (Gabe an Oy 1 — OGG yds (O, — Ogg, DIF 


which leads to 


y[i7] =i V[n] + D- x[77] : 
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(5.34) 
0 
0 
salle pala (535) 
1 
v, [77] 
vl] sees Le 
Vy] 
(5.36) 


Since there are two transfer functions in equation (5.30), we apply the state space 


form twice yielding 


x,[+1]=A,x{/]+ B, ell] 
Y= C, x [+ Del!) 
X,[/ +1] = A, x1] + B, el!) 
Y,olf]=C, xl] + Deli] 


yey) ll yee 


where 


0 
0 0 
A, =A, = 
0 
Nn Ay} 
0 
0 
B, = B, =|. 
] 


C P= Wen Gian of 00... Jy heey 0G) 
C, a [Cf on med 5 eet Sif ce Oe cae Ge a aod.) 


Dy = omane 


D, = fo. 
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Oey) 
(5.38) 
rg) 
(5.40) 


(5.41) 


A, = A, because these matrices are based on the same denominator term D(z7'). 


From the state space model we can now express equation (5.30) in terms of matnx 


variables as 
aCe nc) = +1 
FE ye ath (lata )°B, + Dbl (Elet~A)" B+ Dbl 


(5.42) 


Taking the transpose of expression (5.42) leads to 
yl =(B7 (a -4,7)"C," +, ke, +(B," (a - 4,7)'C," +D, kb []. (5.43) 


Note that D, and D, are scalars. There are clear similarities between the two 


terms in equation (5.43), which is advantageous because we wish to combine the two 


expressions into one state space model. Since A, and A, are the same matrix, the term 


to determine the eigenvalues is equivalent. We combine both terms of equation (5.43) 


resulting in a single state space model for y,[{/], giving 


xl +1) = Axl/] + Bel!] (5.44) 
y,U] = Cxl!] + Del!], (5.45) 


where the matrices in the new state space model are defined as 


Aah =A. 
B = lo oma 
C=B’ =B, 


D=[D, D,],and 
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a ° a 


Looking at the new single set of state space equations (5.44) and (5.45), they 


describe a general stochastic state space model 


alk + 1} = Axtk} + v{k} (5.46) 
z[k] = Cx[k} + w{k]. (5.47) 


An innovation model can be determined from equations (5.46) and (5.47) using a 


Kalman Filter as 


Rk +1] = ARK] + K -ZIK), (5.49) 
Z[k] = 2[k]— CSIk], (5.50) 


where 


z[k] is the observation at time k, 
x[k] is the prediction of state x[k] given all observations up to time k —1, and 


z[k] is the innovation, which is independent on all past observations. 


The Kalman gain matrix K and the covariance matnx P are computed from the 


Riccati equation where 


P = E\X[k)X" [k]} 
with X[k] = x{k] — XK). 


Finally, the covariance of the innovation is given by 
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E(Z{k]’ }=CPCT +R 
where R is the vanance of the white noise process w{[k]. 


For the particular case of equations (5.44) and (5.45), we obtain 


Xl +1) = Axl] + K-e,(1), Goi) 
yi = een (5°52) 


In order to solve the Riccati equation, covariance matrices Q, R and S must be 


determined. Because we assume zero-mean white noise, the covariance matrices are 
identically equal to the correlation matrices. Applying the general state space model 


equations (5.46) and (5.47) with the specific model equations (5.44) and (5.45), we define 
O as 
Q = E{k]-v" (k]}= E{Belle” (1187 J. 


Performing the product of the uncorrelated 2 x N white noise set, we note that 


- ell : 0 il 
elle" (0 -|: net e {|= ° A ade & 


Now substituting back into the above equation and noticing that all the 


expressions in the expectation are deterministic we get 


Q=0°BB’. (5.53) 
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Performing similar operations to R and S$ we end up with 
R=E\ wil] w" [I}=o°DD", (5.54) 
and 


$= E\ vl] w' [t= o°BD". (5.55) 


Note here that the variance o° is the same in equations (5.53) through (5.55) 


because the white noise processes that comprise e[/] are independent, yet have the same 


variance. Recall this variance is determined from the initial ARMA model. 
For the final step in determining the multiresolution innovation process, consider 
the measurement equation (5.52) where we have determined all values in order to 


determine the innovation e,[/] as 
el)= y,[]- Cx{d], 
where it can be shown that the vanance of e,[J] is 
Ce Efe, f= CPC’ +R. (5.56) 


It would be possible to implement this multiresolution innovation process in terms 
of a Kalman Filtering recursion, but notice the state space model representation of 
equations (5.51) and (5.52). This set of equations is also in Type I controllable canonical 
state space form for the Direct Form realization of a transfer function [13]. This means 


that a single transfer function can be determined from the state space model, giving 


y, [1] =[C(zl — A)" K + 1e,[l). >) 
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Finally, it can be shown that if an ARMA model is minimum phase, then all its 
corresponding multiresolution ARMA filters are minimum phase as well. From chapter 
2, recall] that a minimum phase filter will have a causal stable inverse. From this fact, we 
see that the inverse to the multiresolution ARMA model yields our goal of obtaining a 
multiresolution innovation filter bank. 

In a similar way, we can determine a model for the highpass component, as in 
Figure V-1. Consider now the ARMA filter coefficients from the highpass filter channel 
] 


> 


ant 


where G(z_) -| $I. We see that the derivation is quite similar. In fact, the only 


difference is the section where we consider the effect of F(z™') on the Direct Form 
realization of the transfer function in equation (5.42). Going back to equation (5.29) we 


substitute the highpass filter G(z~'), giving 


1-17 Ae) He) Te. 
is |S ee re ° ° D6 
yl ] E g) in (z7") H, es ( ) 


Carrying out the multiplication, substituting in equation (5.16), and grouping 


terms yields the expression 


y= StS n+ S20, (5.59) 
where 

6.0) =Flee)- 2716, (5.60) 
and 


13 


IC, (27) - Oe )| (5.61) 


From this point, we can substitute G,(z') and G,(z~') into equation (5.42), 


which is the state space form for the Direct Form realization of the transfer function in 


equation (5.59). 


B. METHOD 2 SUMMARY 


Method 1 was a method to look at a whitening filter before wavelet 
decomposition. By contrast, Method 2 investigates a method to perform the whitening 
filter after the multiresolution process. The algorithm can be divided into two processes. 
The first one is the standard wavelet decomposition of the orginal signal similar to the 
DWT on the whitened signal for Method 1. The second process in the algorithm for 
method 2 deals with determining the innovation filter bank. A set of ARMA filter 
coefficients 1s necessary for each wavelet decomposition channel, however all] that is 
needed a priori is the ARMA model coefficients of the original signal. That is to say, we 
determine the ARMA filter coefficients based on the original signal and from that model 
we can determine all follow-on multiresolution channels. From this description, we can 
see how the Kalman filtering procedure described in this chapter is applied and a short 
description of the steps for the entire process is shown in Table V-1. 

All the white noise channels are uncorrelated with each other resulting in an 
ability to observe the innovation process with varying degrees of highpass and lowpass 
filtering. By way of analogy, one of the useful properties of the wavelet decomposition is 
denoising. Since most useful signal information occurs in lower frequencies for 
underwater acoustic applications, the wavelet decomposition can leave out some highpass 


information (i.e. details) upon synthesis thereby removing noise from the signal. 
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In our particular application, a similar method of running a signal through a 
constant-Q filter takes place, only now we take each channel and process it through a 
whitening filter. In effect, each channel has less and less high frequency content. It 
would then stand to reason that this multiresolution whitening filter bank might perform 
two desirable functions. First, the filter bank might be able to distinguish low frequency 
narrowband signals well due to its constant narrowing low pass filter, providing high 
resolution in low frequencies. Second, since the channels are uncorrelated, indication of 
a transient signal might be evident in one channel and not in others, depending on its 


frequency content. 


1 


STEP DESCRIPTION EQUATION MATLAB FUNCTION 
Determine ARMA coefficients of original one [, a| =Criste arty) 
1 | signal and variance o° 
Decompose a filter transfer function into SD), 5, 1ke (Cer D]= polyd(num, den) 
Z even and odd parts for Low Pass Filter 
ce) : 
Determine F\(z™') and F.(z7') 5.31, 5.32 ie eee NABE eee 
: model2(b, ,a,) 
Determine the Type I Controllable 565, SoS AB em ay |= 
4 Canonical Form for H,(z™') andH,(z”’) tf 2ss(F,,D) 
Combine even and odd processes to one 5.44, 5.45 bo b rae 
5 State space model model2 (b, ,a,) 
8 
9 
| 













Determine Q, R and S$ Sye)e) 8). Q=0°BB',R=o0°DD' 
S=o’°BD' 
Determine Kalman Gain and State [14] Ox — 
Covariance Matnx from Riccati Equation mydlge(A,G,C,O,R,S) 
Determine new varianceo” ED) o° =CPC' +R 
| 
| | 
Sei 





Determine next multiresolution filter [num, den| = 
transfer function ss2tf (A, K,C,1) 


| 10 Repeat steps ? ~ 12 for HPF G(z") 5 OUe5. al Included in Sler 5 
- 11. | Perform 2 — 13 to the desired number of op ae meth_0102 
Stages 


12 | Run the signal through the multiresolution Figure V-1 meth_O102 
whitening filter bank as in Figure V-1 


Table V-1 Process to Determine the Multiresolution Whitening Filter Bank 
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VI. RESULTS 


Now that we have set a framework for two methods of multiresolution filtering, it 
is necessary to compare the techniques against the benchmark signal described in Chapter 
1. Method 1 can be described as pre-whitening a signal before multiresolution analysis, 
which is the Discrete Haar Wavelet decomposition. By contrast, a signal processed 
through Method 2 uses the same two processes in opposite order. Again, the benchmark 
signal is a fan blowing, which 1s representative of the colored noise environment of the 
heavily — trafficked coastal waters because it contains more of the lower frequencies due 
to a rotating blade through a fluid. 

The test signal used in this section is the first three seconds of the benchmark 
signal unless otherwise specified. The assumption is there are no other transients in this 
regime except for the synthetic signal that starts after the first second. A wide range of 
synthetic transients is tested from narrowband single frequency sinusoids to broadband 
white noise. The object of this section is to compare established detection techniques 
with Methods 1 and 2 using a wide range of transients. 

Since the STFT is the current standard for underwater acoustic analysis, all results 
will begin with a spectrogram followed by the whitening filter process. In this way, the 
viewer Can receive an intuitive “feel” for comparing methods. Introducing the new 
detection methods, as we shall see, results in a higher confidence in positively identifying 
the transient. 

Finally, for ease of comparison, Method 1 and 2 results will be presented 
concurrently at each of the three transient tests. Because of the multiresolution process, 
each method has multiple channels to display. Since not all channels will indicate the 
presence of the signal, all figures left out are assumed to be clutter. For convenience, a 
small diagram of each process will precede every figure with its respective filter path 
highlighted. Knowing which channel is being displayed is important when analyzing the 


effect of the filtering process on the transient signal. 
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A. BENCHMARK TEST 


The method companion will test three different transient types. The first two are 
narrowband transients at opposite sides of the frequency spectrum while the last test will 
be the more general broadband type. We will start with examining the narrowband 


transient signals and make adjustments to Methods 1 and 2 as necessary. 


1. Narrowband Transient: High Frequency Sinusoid 


Note that the maximum testable frequency is half the sampling frequency (just 
over 5 kHz) as the benchmark signal is sampled at 11.025 kHz. Here, we arbitrarily 
choose a 5 kHz sinusoidal transient with a duration of .1 second. The spectrogram in 
Figure VI-1 shows the transient 1s emerging from the background but is barely 
recognizable audibly because the amplitude is small. From this signal representation it 
can be seen that since the benchmark contains little frequency in the region above 2.5 
kHz, a transient signal above this frequency can be easily recognizable, assuming the 
amplitude is high enough. 

Figure VI-2 shows a simple diagram of Method 1 in order to keep track of the 
multiple outputs in the analysis tree followed by Figure VI-3, which shows the innovation 
process on the benchmark signal. Again, in this representation the transient is noticeable 
between the vertical dashed lines. Notice, though, the mean and variance of the noise 
floor compared to Figure VI-5, which is the first stage of the multiresolution process. In 
this procedure, the whitened signal is now processed through a lowpass filter and a 
highpass filter, which are the Haar Wavelet coefficients. 

Looking at the results of the first stage of the Haar Wavelet decomposition shown 
in Figure VI-5, the signal through the lowpass filter provides no useful information, but 
the highpass filter shows the transient very clear and with less variance on the noise floor. 

As a minor digression, it is important to recall that the signal passed through the 
Haar Wavelet filters will result in the Haar Wavelet coefficients, which satisfy Parseval’s 
theorem similar to the properties of the coefficients determined from the Founer 


Transform. In this application, the actual numerical power given for each scale is not 
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important, only their relative values and so all the signals can be normalized. For this 
application, the fundamental properties of the Discrete Wavelet Transform are only 
important in terms of defining the filtering process. To effectively analyze the results, we 
simply view the multiresolution filtering process in terms of how the filters effect the 
signal in the frequency domain. What we have shown in Figures VI-3 and VI-5 is a 


denoising process, which is one of the many uses of wavelet transforms. 
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Figure VI-1 Spectrogram of Transient 1: 5 kHz Sinusoid 
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Figure VI-2 Method 1 Analysis Tree (w[n]) 
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3 Transient 1: Whitening Filter 


Figure VI 
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Figure VI-4 Method I Wavelet Analysis Tree (A, and D., ) 
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Figure VI-5 Transient 1: Method 1, Stage 1 (A, and D,) 
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For the high frequency sinusoid, the transient will not appear in any follow-on 
stages of the lowpass filtering processes (approximations). At this point, we are left with 
a successfully detected signal in the high pass channel of the first stage. For this 
scenario, we can see the transient above the noise floor, but recall this signal can also be 
detected from the spectrogram. Now, suppose we can reduce the amplitude of the 
transient signal to the point where it cannot be easily identified in either the spectrogram 
or the whitening filter. Can the denoising process of Method 1 provide better results? 
The answer is shown in the next set of figures. 

The key to successful denoising on the high frequency end of the spectrum is to 
filter the details (high pass channel) in the same manner as the approximations (low pass 
channel). Knowing that we can remove the signal from the noise on the low pass 
channel, we apply the same wavelet decomposition at every output. This description is 
known as the standard Wavelet Packet framework, which has been applied to a wide 
range of applications such as denoising and data compression. Wavelet packet 
processing allows for a more complex and flexible analysis because both the details 
(highpass channel) and the approximations (lowpass channel) are decomposed, as shown 
in Figure VI-7. In terms of wavelet analysis, the wavelet packet framework allows for a 
wider range of bases from which you can choose the best representation [15]. 

For our application, the flexibility of the wavelet packet allows for greater 
narrowband detection sensitivity in both high and low frequencies. Take the high 
frequency transient for example. Let’s reduce the amplitude of the 5 kHz sinusoid from 
.005 to .003. To quantify this reduction, audibly detecting this transient signal is like 
straining your ears to listen for the softest high frequency sounds during an audiogram 
test. On top of that, we still have the noise of the fan blowing. Figure VI-6 shows the 
spectrogram of the benchmark signal with the barely noticeable transient at 5 kHz, which 
could be classified as noise. Now we add the high pass channel decomposition and find 
that several channels have successfully reduced the noise floor so that the signal can be 
detected. Because of the branching effect of the wavelet packet, a three-stage framework 


will yield 15 plots. In order to reduce the number of plots, only the graphs that detect the 


signal are plotted. Figure VI-8 shows the whitening filter representation of the 
benchmark signal. Notice that the signal is completely imbedded within the noise such 
that it cannot be identified. Performing the first stage of the wavelet decomposition as 
before yields a slight indication of the transient (Figure VI-10), but now we perform the 
filtering process on the details and notice the signal clearly emerges from the noise floor 
(Figure VI-12). Performing the filtering process one more time yields an even better 
result (Figure VI-14). Realistically, the transient might have been positively identified on 
the spectrogram. However with the wavelet denoising process, the added signal detection 
capabilities might add to the measure of confidence. 

At first glance, it would appear that the 5 kHz signal should be detected in the bin 
to the far nght marked DDD, because the bandpass filter in this region contains the 
frequency of the transient. Note in Figures VI-11 to VI-14 that the signal appears in the 
bandpass filter region AD, and AAD,. The reason for this could be due to aliasing that 
comes from the downsampling process is subject for further research. Also, since the 
highpass and lowpass filters of the Haar wavelet do not have a sharp dropoff at the cutoff 


frequency, it is possible to have part of a high frequency transient signal to appear in the 


low frequency channel and vice versa. Both these issues are subject for further research. 


83 





~N 
=) 
> 
S 
Z 
es 
S 
2 
ae 


195 
TIME (seconds) 





Figure VI-6 Spectrogram of Modified Transient 1: Reduced Amplitude 
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Figure VI-7 Wavelet Packet Framework (w[n]) 
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Figure VI-8 Modified Transient 1: Whitening Filter 
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Figure VI-9 Wavelet Packet Framework (A, and D, ) 
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Figure VI-10 Modified Transient 1: Stage 1, Method 1 (A, and D,) 
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Figure VI-11 Wavelet Packet Framework (AD, and DD,) 
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Figure VI-12 Modified Transient 1: Method 1, Stage 2 (AD, and DD,) 
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Figure VI-13 Wavelet Packet Framework ( AAD, and DAD, ) 
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Figure VI-14 Modified Transient 1: Method 1, Stage 3 (AAD, and DAD, ) 
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Turing now to Method 2, recall that the input signal is whitened after the 
multiresolution Wavelet decomposition as described in Figure V-l. Throughout the 
derivation, we used a slightly modified filter bank from the Haar Wavelet coefficients as 
described in equation (5.7). As mentioned before, the only difference is a constant factor 
and since the results are normalized, the shape of the signal is exactly the same. For the 
sake of consistency, the actual Haar Wavelet filter coefficients were used. 

Unfortunately, Method 2 did not detect the narrowband high frequency transient 


signal at amplitude .005. Looking at the lowpass filter, notice that it averages two 


adjacent data points and weights them by a factor of V2. For the case of the 5kHz 
signal, perhaps the transient cannot be seen in the lowpass channel because of its high 
frequency content and cannot be seen in the highpass channel because the transient is 
averaged in with the noise. Nonetheless, to show that the method is still capable of 
detection, it can be shown experimentally that if the amplitude of the transient signal is 
increased to .01, the transient signal emerges from the noise. 

Although the multiresolution innovation process did not perform as well as the 
wavelet packet denoising process at 5 kHz, it is still a viable detection scheme. Indeed, 
we shall see other examples where Method 2 provides the best transient detection 


capabilities. 
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2 Narrowband Transient: Low Frequency Sinusoid 


At high frequencies, the spectrogram works quite well in detecting narrowband 
transients because the stationary colored noise model has very little high frequency 
content. Therefore any added frequency content is easily recognizable. Now consider a 
transient signal buried in the low frequency band between O and 300 Hz. Since the 
colored noise signal has its frequency content in this range, a signal cannot be 
distinguished from the background noise. Since both Methods 1 and 2 are based on a 
signal’s statistics in the ttme domain and not on the Founer Transform in the frequency 
domain, they should detect low frequency transients equally as well. At each stage of the 
wavelet decomposition, the low frequency channel continues to be passed through a 
lowpass filter, effectively cutting the frequency content of the remaining signal in half. 
Continuing this process provides infinite resolution in the low frequency channels, which 
again is one of the nice properties of the Discrete Wavelet Transform. For this reason, 
Method 2 should perform at least as well as Method 1 in the high frequency narrowband 
case. 

We choose an arbitrary sinusoidal signal at 50 Hz and conduct a similar analysis. 
Figure VI-15 shows the spectrogram of the benchmark signal along with the imbedded 50 
Hz transient signal. It is hard to tell whether or not the transient can be detected in this 
region because it could easily be identified as a spurious point. The transient clearly 
cannot be detected with just the whitening filter (Figure VI-17) and not even audibly. 
Starting with Method 1, we see the signal emerging on the low pass channel in Figure VI- 
19. Because the frequency is so low, the signal becomes more pronounced with each 


successive lowpass filter as shown in Figures VI-20 through 23. 
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Figure VI-15 Spectrogram of Transient 2: 50 Hz Sinusoid 
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Figure VI-17 Transient 2: Whitening Filter 
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Figure VI-18 Wavelet Packet Framework (A, and D,) 
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Figure VI-19 Transient 2: Method 1, Stage 1 (A, and D,) 


95 


DA, 


Figure VI-20 Wavelet Packet Framework (AA, and DA,) 
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Figure VI-21 Transient 2: Method 1, Stage 2 (AA, and DA,) 
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Figure VI-22 Wavelet Packet Framework ( AAA, and DAA, ) 
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Figure VI-23 Transient 2: Method 1, Stage 3 (AAA, and DAA, ) 


Method 2 results for the low frequency signal are just as revealing. Unlike the 
Wavelet Transform on the white noise process, the colored noise signal does not have 
equal frequency content at all frequencies (except for the transient) and therefore the 
narrowband transient will not fall on one side of the highpass/lowpass filtering process. 
Indeed, the signal may show up in as many as all of the channels as 1s the case here, 
which in and of itself might provide added redundancy since the channels are 
independent of each other. For the lowpass-filtering channel, the benchmark signal is 
downsampled and averaged. As it averages adjacent points, the process essentially 
reduces the noise level while maintaining a strong indication of the presence of a signal. 
For the highpass filter, the difference between adjacent points is computed and the 
combination of the two processes form a basis set for the onginal signal. Again, this 
highpass process might prove valuable since it is independent of the lowpass process. 
Figure VI-25 shows the signal clearly emerging in both channels. With each successive 
multiresolution innovation process, the noise floor is reduced and the signal becomes 


more pronounced (Figures VI-26 to 29). 
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Figure VI-24 Multiresolution Innovation (A, and D,) 
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Figure VI-25 Transient 2: Method 2, Stage 1 (A, and D, ) 
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Figure VI-26 Multiresolution Innovation Stage 2 (AA, and DA,) 
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Figure VI-27 Transient 2: Method 2, Stage 2 (AA, and DA,) 
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Figure VI-28 Multiresolution Innovation Stage 3 (AAA, and DAA,) 
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Figure VJ-29 Transient 2: Method 2, Stage 3 (AAA, and DAA,) 
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3. Broadband Transient: Elvis 


We have studied the effects of Methods 1 and 2 on narrowband transients. In the 
case of underwater acoustics, it is quite possible that the signal we wish to detect is a 
single frequency like an active sonar pulse or a unique single frequency tonal. On the 
other hand, a broadband signal such as the splash of a mine hitting the water or the clang 
of a metallic object 1s also possible. 

From the narrowband transient, we now consider a more general case where the 
expected signal has a wide frequency band. For this test, the frequency content must be 
completely arbitrary, so for no apparent scientific reason the arbitrary colored noise 
signal will be Elvis singing that great Christmas classic, “Blue Christmas’. Besides, it is 
just as likely to find Elvis in the deep blue ocean as it is to find him anywhere else. 

In order to see the frequency content of the narrowband transient signal, a 
spectrogram of the transient alone is provided in Figure VI-30(a). Notice that Elvis can 
really hit those low notes. The strength of the frequency content is large in the range 
between 0 and | kHz, then tapers off as it approaches 5 kHz. Adding Elvis to the colored 
noise model as in Figure VI-30(b), notice that the signal can barely be seen, and that only 
because we have other spectrograms of the same data with which to compare. Without 
the redundant benchmark data, “Blue Christmas” is reduced to mere noise. We perform 
the innovation process on the benchmark signal in Figure VI-32 and, no doubt about it, 
Elvis is clearly alive. Performing the denoising on the innovation as in Method 1 
confirms this to be true, however the signal is less pronounced. Also notice that the 
signal shows up — in part — on both highpass and lowpass channels at the same stage. 
This indicates that the broad range of frequency content of Elvis’ golden voice is divided 
into highpass and lowpass regions. The problem with Method 1 here is that along with 
the reduction of noise with each filter stage comes the reduction of the signal power as 
well depending on its frequency content. In other words, the success of Method | 1s now 
dependent on frequency content. 

The signal had the same effect when processed through Method 2, only worse. 


Only the first stage provided any indication of Elvis’ location as in Figure VI-44. Here 
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again, the broadband signal could have been averaged in with the noise so that the 


innovation process could not distinguish Elvis from a fan. 
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Figure VI-30 Spectrogram of Elvis Without Benchmark Colored Noise 
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Figure VI-31 Spectrogram of Transient 3: Elvis 
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Figure VI-33 Elvis Through a Whitening Filter 
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Figure VI-34 Wavelet Packet Framework (A, and D,) 
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Figure VI-35 In Search of Elvis: Method 1, Stage 1 (A, and D,) 
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Figure VI-36 Wavelet Packet Framework (AA, and DA, ) 
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Figure VI-37 In Search of Elvis: Method 1, Stage 2 (AA, and DA,) 
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Figure VI-38 Wavelet Packet Framework (ADA, and DDA, ) 
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Figure VI-39 In Search of Elvis: Method 1, Stage 3 (ADA, and DDA, ) 
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Figure VI-40 Wavelet Packet Framework (AD, and DD,) 
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Figure VI-41 In Search of Elvis: Method 1, Stage 3 (AD, and DD, ) 
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Figure VI-42 Wavelet Packet Framework (ADD, and DDD,) 
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Figure VI-43 In Search of Elvis: Method 1, Stage 3 (ADD, and DDD,) 
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Figure VI-44 Multiresolution Innovation Stage 1 (A, and D,) 
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Figure VI-45 In Search of Elvis: Method 2, Stage 1( A, and D,) 
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B. SUMMARY OF OBSERVATIONS 


We have presented a representative set of transient signals tested on two detection 
schemes. This section aims to provide a broad brushstroke of the effectiveness of the 
Methods as they are compared to the spectrogram and the whitening filter. To 
accomplish this comparison, a measure of successfulness must be defined. Since the 
location and duration of the transient signal is known, an estimate of the signal power and 
the noise power can be calculated as the Signal-to-Noise ratio. Since our signal model 
has additive noise, the signal power is simply the power above the noise threshold 
divided by the noise power where both power terms can be calculated from the whitened 
signal. For Methods 1 and 2, whose multiple channels will not always provide the best 
transient detection, only the channel with the maximum SNR will be taken. The 
assumption here is that there is a way of always determining the best channel within a 
given method scheme. 

For a signal whose power is twice the noise power, the SNR will be 3 dB. This 
value will be our threshold detection mark where any SNR above this value is assumed to 
have positively detected the transient. Any value below this mark is assumed to not be 
able to identify the transient even though a visual representation might indicate 
otherwise. 

As in all experiments, the SNR measure is only an approximation. In our case, 
the same exact transient signals within the colored noise model is used to compute the 
SNR and so the results among the different methods may be compared to one another. It 
is fair to assume, however, that the results from this experiment are indicative of many 
types of transient signals hidden within an arbitrary colored noise scheme. The first 
experiment deals with broadband signal detection followed by a test on narrowband 


signal detection. 


lekO 


le Benchmark Comparison 


The first experiment begins with Elvis. Taking the signal, we increase the 
normalized volume from 0 to 1 and determine its SNR from the whitening filter alone, 
Method 1 and Method 2. Figure VI-45 shows the results. Here, it can be seen that 
Method 1 reaches the 3-dB threshold first, followed closely by Method 2 and the 
benchmark whitening filter. As Elvis gets louder, though, Method 1 reaches a higher 
maximum SNR value of 25 dB while the other two methods are just about identical as 
they reach a steady state value of about 20 dB. The results of the experiments show that 
the noise level does indeed reduce faster than the signal and so a better detection method 
was found even for arbitrary broadband signals. On the other hand, comparing any of 
these results with the threshold reveals that Method 1 might reach detection status before 
the other methods, but not by much. What’s more, the spectrogram might be able to pick 
up the signal before any of the three signals, but again there must be some frequency 
content outside of the frequencies defined by the colored noise. 

Another broadband test was conducted where the frequency content expanded the 
entire spectrum. A similar test was conducted where the broadband signal was white 
noise and the results are found on Figure VI-46. Here again, Method 1 broke out with the 
highest maximum SNR at 40 dB. Since all three methods reached the 3 dB threshold at 
about the same amplitude, all three methods perform about the same. Also, we are 
guaranteed to have the spectrogram detect the transient at some point close to the other 


methods and so any of these processes are equally effective. 


11) 





Threshold 


- Whitening Filter 
Method 1 
Method 2 


jt =n —L 


02 03 04 #05 06 O07 08 09 
AMPLITUDE OF ELVIS 





Figure VI-46 SNR Comparison on Colored Noise Transient 








Threshold 


| - W hitening Filter 
Method 1 
Method 2 — 


02 03 04 05 06 07 08 09 
AMPLITUDE OF BROADBAND SIGNAL 








The final comparison is with the narrowband transient. For this test we vary both 
frequency and amplitude of a standard .1 second transient and process this signal through 
the whitening filter, Method 1 and Method 2. The resulting 3-dimensional image reveals 
some interesting peaks and valleys. Figure VI-47 shows the result from the whitening 
filter. Comparing this plot with the results from Method 1 (Figure VI-48) we see nght 
away that it is the same basic shape, only the SNR has increased by about the same 
amount over the entire frequency-amplitude plane. Now recall that the threshold is 3 dB 
and notice how well the wavelet denoising process worked. In the same manner, Method 
2 produced very good results in the low frequency region. 

We wish to further quantify the results of the comparison test. Taking the 3- 
dimensional plot, we slice the graph along the frequency-amplitude plane at a constant 
level of 3 dB. Looking down on the slice, we can see clearly the regions where the SNR 
is above 3 dB, which is the detection region. Figure VI-49 shows the modest narrowband 
transient region for the whitening process. With the exception of a small strip, the lower 
frequency SNRs remain below the threshold, especially between 0 and 300 Hz where 
there is no SNR above 3dB. Now Comparing Method 1 to the whitening filter detection 
region (Figure VI-50) we notice nght away the larger areas in which the SNR is above 
the threshold. To quantify, Table VI-1 shows that only about 30% of the entire 
frequency-amplitude plane is above the threshold. With the Haar Wavelet denoising 
method, the percent detection over the entire region increased to 74%. In the low 
frequency region between O and 2 kHz, the detection region increased from 29.7% to 
78.3%. The results from Method 2 as seen in Figure VI-51 show an equally strong 
performance in the low frequency regions and as mentioned earlier, it has limitations in 
the high frequency regime. 

These results seem very good compared to the whitening filter, now compare the 
methods against the spectrogram. The spectrogram can detect most signals above an 
amplitude of .1 and frequencies above 300 Hz. Recall that the whitening filter did not 
have any of its detection region below 300 Hz so with these two methods combined, low 


frequency transients would go unnoticed. Looking back at Figures VI-50, VI-51 and 


Table VI-1, we notice that the detection region above the 3 dB mark is now almost 50% 


of the total area in this region. 
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PERCENT DETECTION IN FREQUENCY RANGE 


FILTERING 
TECHNIQUE 0 - 5 (kHz) | 0-2 (kHz) | 0- 300 (Hz) 


MULTIRESOLUTION 
INNOVATIONS 





Table VI-1 Narrowband Transient Detection Capability Comparison Over a 


Threshold of 3 dB 
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Figure VI-48 Narrowband Transient SNR For Whitening Filter Process 
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Figure VI-49 Narrowband Transient SNR for Haar Wavelet Denoising 
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Figure VI-50 Narrowband Transient Detection Region For Whitening Filter 
(Threshold = 3 dB) 
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Figure VI-51 Narrowband Transient Detection Region For Haar Wavelet Denoising 


(Threshold = 3 dB) 
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Figure VI-52 Narrowband Transient Detection Region For Whitening Filter 
(Threshold = 3 dB) 
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Figure VI-53 Narrowband Transient Detection Region For Multiresolution 


Innovation Process (Threshold = 3 dB) 
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Zz, Strengths 


The most useful property of the wavelet denoising process and the multiresolution 
innovation process are their abilities to detect narrowband signals in the low frequency 
regime, particularly in between O and 300 Hz. This property is highly desirable since it 
has been shown that neither the spectrogram nor the whitening filter is effective in this 
frequency region. | 

The Haar Wavelet denoising routine works well at high frequencies too because 
of its wavelet packet framework. Here, the high frequency details are filtered too, 
creating greater flexibility by isolating more frequencies. 

Since the multiple-channel outputs of both Method 1 and Method 2 are 
independent of each other, it is likely that more than one channel will indicate the 
presence of a transient signal. This property provides a measure of redundancy in 
positively identifying a transient signal. 

Both methods are simple to compute numerically; the DWT is computed by a 
simple filtering process as the number of computations 1s comparable to current methods 
that use the FFT. Since the size of the data set is reduced by a factor of 2 at each stage, 
the number of computations is reduced by the same amount. This property is desirable 
because we know that any follow-on stages will incur less and less computations. 

Both methods lend themselves to an efficient modular processing routine with a 
user-friendly GUI display. Take for example the wavelet packet framework. The user 
might designate the number of stages to explore and at any time, can add another stage 
without having to go back and compute the entire wavelet tree. The same is true for the 
multiresolution whitening process because, like the wavelet decomposition, every 
successive ARMA filter set is only based on the previous stage and does not need to be 


computed from the very first stage again. 





S Weaknesses 


Both wavelet denoising and multiresolution innovation procedures attempt to 
reduce the noise variance by a series of filtering processes. For a signal that is isolated to 
a small range of frequencies, both methods will perform quite well since the filtering 
process has the ability to cut off a sizeable amount of extraneous frequency content. By 
contrast, a broadband signal will have part of its frequency content cut off at each 
filtering stage. The effect will be a reduction of signal power along with a reduction in 
noise power. Even with the signal power reduction, the results of the previous section 
show that Methods 1 and 2 were able to reduce the noise floor as the amplitude of the 
transient signal increased. Comparing these methods to the spectrogram, however, we 
see that the transient will appear as a vertical line at amplitude less than .05, which is 
about when the signal emerges from the noise in the multiresolution filtering process. 
The advantage of the spectrogram for broadband signals is its color representation, which 
essentially adds an extra dimension to detecting the transient. The only exception to this 
observation 1s when the transient contains all or part of the same frequencies as the 
background noise, in which case the signal is camouflaged in the spectrogram. 

Inherent in the wavelet decomposition process is a reduction of resolution in the 
time domain. As described in Chapter 3, to perform a DWT on a signal it is 
downsampled as part of the process. By this action, the signal in the time domain loses 
some detail at each stage. Figure VI-53 shows the entire 6 seconds of the benchmark. 
Recall the peaks located in the latter half of the benchmark signal are a set of 11 
broadband transients. Notice the cluster of pulses just before the fourth second in the 
graph of the highest sampling rate. Here, the three spikes are easily recognizable. Now 
dropping down through the other three graphs where the sampling frequency is reduced, 
notice that the three spikes are slowly combining to form a single transient. If the final 
stages of method 1 or 2 were the only channels to pick up this set of transients, there 


would be no way to differentiate between them. 
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Figure VI-54 Reduction of Signal Resolution Due to Downsampling 


VII. CONCLUSIONS 


This thesis introduced new transient detection techniques designed to provide 
sonar operators an enhanced capability to detect enemy acoustic signals. The methods 
presented in this thesis are based on the Autoregressive-Moving Average and Discrete 
Wavelet Transform. The difference is the order in which these two operations are 
performed. Method 1 first whitens the input signal using an ARMA model before 
denoising the whitened signal with a wavelet decomposition operation. Method 2 first 
applies wavelet decomposition before applying a bank of whitening filters to each 


wavelet channel. 


A. USW LITTORAL WARFARE APPLICATIONS 


The experimental results reveal a significant increase in signal detection using 
either proposed method when the transient is narrowband. This ability to localize 
frequencies might be useful for applications where the frequency of a potential transient 
is known a priori, such as for example in active sonar operations where the return signal 
is known to be within a Doppler shift of the transmitted frequency. Wavelet denoising 
has the ability to narrow its filtering process such that a specific channel can be observed, 
rejecting all other frequency content. Therefore, this process has the potential to increase 
the effective range of the sonar without having to increase the signal power. 

Another application to narrowband signal detection might be for passive sonar 
operations where one attempts to uncover a known frequency within noise of similar 
frequency content. In low SNR, we have shown that the spectrogram does not provide 
useful capabilities in this regime and so we rely on time-domain methods. Examples of 
such signals might be specific tonals or blade rates produced by enemy platforms hidden 
within the noise of background shipping. 

In spite of the fact that a computer can successfully detect a transient signal over 


some given threshold, the sonar operator is still a vital member in the detection phase of 


USW. Indeed, it is not unusual for the initial detection of an enemy submarine to come 
from a source other than an automatic alert function. Compared to the other warfare 
areas, human intervention 1s more vital in the detection phase because the pace is 
generally slower. With this in mind, an operator must have multiple displays and sensors 
at his disposal. It is believed that the methods tested in this thesis will not only increase 
detection capabilities, but will also provide enhanced display features. 

The main idea behind multiple display features at an operator’s disposal is to 
correlate data. Consider a Graphical User Interface (GUI) system where an operator is 
shown a wavelet packet framework diagram as in Figure VI-7 where the boxes represent 
the wavelet decomposition channels. The display feature could also link with existing 
features in the sonar suite to further enhance transient detection. 

Another way to enhance operator displays would be to include various intensity 
levels that somehow correspond to the SNR. Examples include full color displays or 


alarms with a range of volumes relative to the strength of a transient signal. 
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B. SUGGESTIONS FOR FUTURE RESEARCH 


As the wavelet decomposition can theoretically be carried on to infinite detail, so 
can the ideas of improving technology. Although this section is certainly not all- 
inclusive, several suggestions are provided for future research based on the results from 


Methods 1 and 2. 


1. Method 1 Improvements 


The DWT denoising process of Method 1 uses the Haar Wavelet, which is the 
most basic wavelet basis function that uses only a first order polynomial for the highpass 
and lowpass filter coefficients. As one can see from the frequency response of the filter 
coefficients (Figure ITI-9), there is a significant overlap between the two filters, which 
results in a significant amount of high frequency information contained in the low 
frequency channel and vice versa. Further research can be done to study other wavelet 
basis functions that have a sharper frequency response. To see a preliminary example of 
denoising using a different wavelet basis function, refer to Appendix B. 

Because we know that the wavelet decomposition produces independent channels, 
another area for follow-on research is to define a way to correlate data at each channel, as 


described in Appendix B. 


2: Method 2 Improvements 


Recall the successful improvements the dyadic filtering had on the narrowband 
high frequency transients. It is believed that the same success will occur for Method 2 if 
the wavelet decomposition were carried out on the highpass channels (details) as well. 
This would require further investigation on defining the whitening filter coefficients for 
these channels. In effect, Method 2 would then perform the full wavelet packet 
decomposition on the input signal and then run it through the modified innovation filter 


bank. Like Method 1, the modified Method 2 would result in 15 channels for the three- 


stage wavelet decomposition. Again, this structure lends itself to a useful GUI display in 
which access to a large number of signal data is relatively simple. 

In addition to investigating the wavelet packet structure for the DWT-Innovation 
process, the same argument can be made about possible improvements using another 
wavelet basis set such as the Daubechies family. This process is a bit more complicated 
than changing the wavelet set in Method | because the innovation filter bank of Method 2 


is dependent on the wavelet filter coefficients. 


bp General Improvements 


The benchmark signal was a colored noise model with synthetic signals 
representative of coastal water acoustics. The main thrust for future study should be to 
use actual acoustic data in this regime. In this way, a more realistic conclusion can be 
drawn for this particular application. 

Results showed that the denoising process performed well for both Methods 1 and 
2. Perhaps the two methods could be combined such that we take the whitened signals 
from Method 2 and pass them through the DWT denoising process of Method 1. 

Looking at the benchmark signal, note that the three-second window performed 
well for the .1 second transient signal. However, this window size might not work well 
for signals of different length. If we wish to detect a range of transient signal durations, 
say as long as 6 seconds and as short as .01 second, further research must be done to 
investigate the effects of different window sizes with varying transient lengths. 

Finally, the last suggestion for future research would be to investigate methods for 
applying a neural network scheme to classify patterns of signals among noise in order to 
detect transients at a lower SNR. If there is a way to determine features that positively 
identify transient signals, then for USW applications it can be another valid tool to make 


the detection process more robust. 
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APPENDIX A. HAAR WAVELET PROPERTIES 


A. NECESSARY TIME DOMAIN CONDITIONS 


This appendix shows that the Haar wavelet system meets all 8 necessary 
conditions for a valid wavelet system [9]. Starting with the wavelet decomposition tree 
in Figure III-6 and the MRA equation, we wish to determine the Haar coefficients of the 
PR filter bank, the scaling function and the shape of the Mother Wavelet. Based on the 
set of eight necessary conditions, we will show that the Haar system of equations is 
indeed a valid wavelet system. 

Recall that all parameters that define a wavelet system can be determined from 
scaling function. If we can somehow develop a parameterizing scheme to define the 


scaling coefficients h(n) , we could determine any wavelet system. To this end, we begin 


by choosing an arbitrary filter size N. The filter size for the scaling coefficients 


h(n) must be even and must satisfy both the linear constraint 


Sh(n) = V2, (A) 


Noe 
and the = bilinear constraints which leads to 


— 


S|a(n)! =1. (A.2) 


It is important to note that the filter (7) is an FIR filter, which also implies 


compact support. This wavelet property is highly desirable as it aids in the time- 


localizing ability of the DWT and it also reduces the number of computations. 


N 
It can be shown that there are om 1 degrees of freedom (parameters) in choosing 


— 


h(n) that will guarantee the existence of a valid wavelet system [9]. Consider the 


simplest case where N =2. This filter length means that there are zero degrees of 


freedom and by equations (A.1) and (A.2), we get 


130 


h(0) + A(l) = V2, (A.3) 
and 
h?(0)+h?QM=1. (A.4) 


Although the second equation is nonlinear, it can be seen by inspection that the 


unique solution to this set of equations is: 


h(n) = {h(0), h()}= 35 a al (A.5) 


These are the Haar scaling function coefficients, which were developed in 1910, 
separate from studies in modern wavelet analysis. In 1992, Ingrid Daubechies showed 
that Haar’s work on this set of orthonormal basis functions are in fact the simplest set in a 
larger family of wavelet systems. From equation (3.19), we can determine the wavelet 


coefficients as 


h,(n) = {h, (0),fy ()}= 5 7 t (A.6) 


With just the coefficients to the Haar scaling and Wavelet functions, it can be 


shown that this wavelet system meets eight necessary conditions [9]. 
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ie Theorem 1: 


If p(t)€ L’ is a solution to the MRA equation (3.14), and if | p@)at #0, then 


Sih(n) = /2. (A.7) 
For the Haar Wavelet the above equation leads to: 


h(0)+h(1) = + 


l 
3 1B 


o|- 


2. Theorem 2: 


If o(t)is an L’ solution to the MRA equation (3.14) with | o@at #1, then 
Yer-D=SeM=1, (A.8) 


wih O(7 + 27k) #0 for some k, then 


SY) h(2n) = »y h(2n +1), (A.9) 


where (3.30) may have to be a distributional sum. Conversely, if (3.31) is 


satisfied, then (3.30) is true. 
For the Haar Wavelet the above equation leads to: 


h(O)+ h(2) = h(Q) + h(3) 


3. Theorem 3: 


If g(t)is an L* AL’ solution to the MRA equation (3.14) and if integer translates 


of Y(t) are orthogonal as defined by 


E k=0 
[omg —k)dt = Eo(k) = | | (A.10) 
QO otherwise 
then 
1 k=0 
Y A(n)h(n - 2k) = 6, = . (A.11) 
7 QO otherwise 


For the Haar Wavelet the above equation leads to: 


h(O)h(-2k) + A(I)A( — 2k) 
h(0)h(0) + A(1)A(L) = ; + =1 for k=0 


h(0)-0+h(1)-0=0 fork<O and k>O. 


Corollary 1: 


Under the assumptions of Theorem 3, the norm of h(n) ts automatically unity, 


L.@., 


S|h(n)| = 1. (A.12) 


n 
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For the Haar Wavelet, corollary 1 leads to: 


h*(0) +h (I) => 45 =1. 


Corollary 2: 


Under the assumptions of theorem 3, 
Yh(2n) = } h(2n+1) = — 
V2 


For the Haar Wavelet, corollary 2 leads to: 


h(O) + h(2) = h() + h(3) = as 


5 


4. Theorem 4: 


If g(t)has compact support on OStSN-1 and if g(t—-k)are linearly 


independent, then h(n) also has compact support over OSn<SN-1: 


h(n) =Oforn<0O and n>N-1l. (A.13) 


For the Haar Wavelet the above equation leads to: 


Haar Wavelet has compact support as /i(72) = 0 


fom 0 andy a1: 
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B. NECESSARY FREQUENCY DOMAIN CONDITIONS 


The last four theorems are necessary conditions in the frequency domain. To 
begin the analysis in the frequency domain, we look at the frequency response of the 
scaling coefficients h(n). The Haar Scaling filter coefficients in the z- domain form the 


equation 


hy(zt)=—=(1+ 27), (A.14) 


S|- 


Substituting e’” for z we get the frequency response of the FIR filter 


H,(w) Se iesee (A.15) 


15 


Likewise, the frequency response of the Haar Wavelet filter coefficients is 


FAC) —— (leo (A.16) 


Wp 


It is important to note here that the frequency response of the scaling coefficients 
h(n) (also known as h,(n)) is that of a lowpass filter with its cutoff frequency at the 
Nyquist frequency rate (H(7z7) =0). It can also be shown that the frequency response of 
the wavelet coefficients h,(m) is that of a highpass filter with its maximum at 
H(z)= V2 and its zero magnitude at H(0). This is the property that allows the input 


signal to be filtered into a “constant-Q” filter bank resulting 1n equal highpass and 


lowpass portions preceding every downsample as described in Figure III-6. Figure III-8 


shows the frequency response of the PR filter bank used for the Haar Wavelet system. 


1. Theorem 3: 
If g(t)is a L’ solution of the MRA equation (3.14), then the following equivalent 


conditions must be true: 


h(n) = H(0) =~V2. (A.17) 
2 


For the Haar Wavelet the above equation leads to: 


H(0O)= os (1 - 1) = AD , therefore 


aD 
h(0) + h(l) = H(0) = +2. 


fae Theorem 6: 


For h(n)é L’, then 


S'A(2n) = SY A(2n +l)if and only if H(7r) =0. (A.18) 


For the Haar Wavelet the above equation leads to: 


ING = ani —1)=0, therefore 


iD 


(0) + h(2) = hl) + h3) = 


‘| 
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3: Theorem 7: 


If p(t) is a solution to the MRA equation (3.14) in L> AL’ and ®(w) is a 


solution of 
@(w) = z aS 3} (A.19) 


and after iteration becomes 


= i W 
aD) = — H| — |'®@(0), A.20 
(W) | (s} (0) (A.20) 
then 


| g(t)p(t — k)dt = 5(k) if and only if ¥'\e(w+ 2al)° =1. Ga 


This theorem 1s the frequency domain equivalent to the time domain definition of 
the orthogonality of the scaling function [8]. Equation (3.42) shows that for a function 


that is narrow in time 1s spread out in frequency. 


For the Haar Wavelet the above equation leads to: 


@(qw) is the Fourier Transform of the scaling function g(t). Finding ®(q@) is an 


iterative process and equation (3.40) is often called the “Cascade Theorem”. By 


assuming an initial ®(Q), the use of a computer to iterate the process can show that 
@®(qw) converges, its inverse Founer Transform is the scaling function and equation 


(3.42) is satisfied. As stated before, once the scaling function is determined, then by 


equation (3.18) one can determine the shape of the Mother Wavelet. The shape of the 


Haar Scaling and corresponding Wavelet function is shown in Figure III-7. 


4. Theorem 8: 


For any h(nje Li, 


»> h(n)h(n —- 2k) = 6(k)  ifand only if IH (a)" ~ IH (w ~ m\" = 2. (3.47) 


For the Haar Wavelet the above equation leads to: 


Starting with H (@) = hh +e /® I. we take the absolute value as multiplying by 


/2 


its complex conjugate giving, 


| eae Cee oe Gree ae 


1 to -je ° ee — a de J + — 
V2 2 V2 V2 ag Dd 
Knowing that the complex exponential form of the cosine function is 


Cos(x) = (e® + oe) 


Wie 


we substitute this into the above expression giving 


lH (w)| = leona) 


Now taking H(w+7) = ( + e storm) ) and multiplying by its complex 


S]|- 


conjugate gives, 
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— 
.. 
Se 
® 
tne, 
a 
-$ 
® 
t 


! a =- 4 e 
2 


Bip late 
V2 2 V2 2 


substitute the cosine function in for its complex exponential form 
The 


=-1., 


Again we 
except now it is shifted by 180 degrees because of Euler’s Identity e— 


resulting expression yields 


IH(@ +7)” =1-cos(w). 


Now we add both expressions together resulting in the final expression 


IH (w)| +|H(w+2)° =1+cos(w) +1-cos(@) = 2 


I? 


THIS PAGE INTENTIONALLY LEFT BLANK 
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APPENDIX B. FURTHER RESEARCH 


The denoising procedure used in Method 1 uses the Haar Wavelet system, which 
is the most basic wavelet set. Looking at the wavelet decomposition as an inner product 
between the basis set and the signal, we see that the Haar wavelet (Figure III-8) works 
well with functions that are piecewise continuous, but not smooth as in the case of most 
signals. Therefore it would be better to find a basis set where the Mother Wavelet has a 
shape similar to that of the signal, such as a wavelet from the Daubechies family [15] 
shown in Figure B-1. By denoising using a smoother, compact wavelet we still get the 
highpass and lowpass effect, but now the wavelet coefficients (output of the filter) might 
be larger due to better correlation between the signal and wavelet. Consider the Haar 
Wavelet in the frequency domain shown in Figure III-9. It has been shown that the Haar 
filter coefficients are only first-order polynomials. These basic filter coefficients imply 
that the cutoff frequencies of the highpass and lowpass channels overlap significantly. 
Therefore, the DWT denoising process might perform better with a filter pair whose 
frequency response is sharper at the cutoff frequency like Daubechies 10 (db10) as shown 
in figure B-2. Indeed, upon one simple test using the same 5 kHz transient signal, the 
db10 wavelet system was able to reduce the noise threshold lower than the Haar wavelet, 


as shown in Figures B-3 through B-6. 


14} 








Figure B-1 Daubechies 10 (db10) Wavelet 
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Figure B-3 Wavelet Packet Framework (AAD, and DAD, ) 
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Figure B-4 Denoising 5 kHz Signal Using Haar Wavelet System (AAD, and DAD, ) 
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DA, : DD, 


Figure B-5 Wavelet Packet Framework (AAD, and DAD, ) 
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Figure B-6 Denoising 5 kHz Signal Using db10 Wavelet System (AAD, and DAD, ) 
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Another modification to Method | that might be considered for further research is 
to define a Figure of Merit (FOM), which would allow for comparison of the various 
wavelet channels. To this end, we investigate the statistics of white noise. By definition, 
a zero mean white noise process x(t)is defined by its first moment (or mean) described 


as 
Ejx(t)}=0 (B.1) 

and second moment (or variance) described as 
E{x(1)x(t -7)}=0,,78(r). (B.2) 


In general, we wish to come up with an operation that can describe the correlation 
between N white noise processes, where N is the number of channels in the wavelet 
decomposition and the resulting correlation matnx R,,. is of dimension NxN. In 


addition, this operation must provide visual clarity in distinguishing a transient signal 


from background noise. With this in mind, a figure of merit is developed as 





min eig(R,, ) 


FOM =1- (B.3) 


max eig(R,., ) . 


The ratio of the minimum eigenvalue to the maximum eigenvalue will always be 
a number between zero and unity. If the minimum eigenvalue is equal to the maximum, 
indicating a set of completely uncorrelated signals with normalized power, then the ratio 
would equal 1 and the whole FOM would equal zero. On the other end, if there is some 
correlation between any two signals (i.e. a transient detected in at least one channel of 
the multiresolution filter bank), the ratio in the expression (B.3) will be less than unity 


and the entire expression will be greater than zero. This non-zero value will show up as a 


145 


“spike” when the FOM 1s plotted as a function of ttme. The square root simply smoothes 
out the FOM so that an indication of a transient signal is easier to see on a visual display. 

There is one final note about this method. Because of the downsampling 
operation after the filter bank, the length of data will decrease by a factor of two at each 
stage. In order to make every vector the same length, the signal is upsampled and run 
through an averaging filter. The effect of this process “stretches” the signal to the desired 
number of samples, and makes every other sample the average of its two neighbors. 
Performing this operation does not affect the statistics of the signal, nor does it wash out 
any transient information. The other option to perform the correlation operation is to 
downsample every detail to the size of the finest approximation. This too would be a 
valid process, but with downsampling comes the loss of information. The other reason to 
upsample and average is that the correlation function determined by a computer is an 
averaging operation. The upsample-and-average method provides more data points to 
average and therefore a better correlation representation. After calculating the details and 
approximations, the data 1s run through an upsample/averaging process to achieve the 
proper data length (indicated as primed details and approximations in Figure B-7). 

The re-sizing of the detail and approximation at the finest scale are then combined 


together as a set of column vectors as shown in Figure B-8 and in the equation 


(lead 
Dee iis dial (B.4) 
as. 
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Figure B-8 Modified Method 1 Display 
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For this procedure, three stages of multiresolution are again used. Finally, the 


correlation function is computed as 


R. =—w -w, (B.5) 


and the figure of merit is determined as 





min eg 5S a (t) 


A plot of the Figure of merit reveals values going from zero to unity as correlation 


between any two inputs increases. A “spike” on the FOM plot reveals a transient signal. 
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APPENDIX C. MATLAB PROGRAMS 
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CODE COMMON TO BOTH METHODS 


| ARMA model 


EUNC Gren (S,el=ecrrsunuar (v7); 
% 

% |b, al=eCristtyasiy 

% 

% vy = data (a COLUMN vector) 


PV=CYistrecosaiyy ie 


$ Levinson recursion 
a=l1; b=l; sSig2=ry(1); sig2p=sig2; 


p=0; £act=0; 


%$ while (p<i0); 

wollte (p<10)4(taect<0.98); 

Bia es +2 Jie 

Pair Timuie (a)v- Dperime=4 = 11 pud Gor 
g=D/Sig2p; gprime=Dprime/sig2; 
ap=(4a;0]—g" (0; flipud(b) i; b=(5; 0)=gprime* (0: tlipud(a) |- 
a=ap; 

fact=l-g*gprime; 

Sige- Pact sige. 

Sige p=face™Sig2p, 

p=p+1; 

end 


b=sqrt(sig2p); 


2 Autocorrelation function 


ED aG elon Steve—tene sesh aloparg (ee) 

S$Function will calculate the cross-correlation function. 
Tt 1s faster than MATLAB function XCORR 

6 Gh xy =Cr Seiecort is 4) 


m=max(({length(x), length(y)]); 
EM=fEt (x) 2 Ie 

EV=ELE(yY, 2m 
SxV=£xX.. * COnnuer 7); 
Rxy=real(ifft(Sxy) )/m; 


< 


send SuCr Use mecorr 
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3% Benchmark Test 


$Name: John D. Stevens 
SDatcrmeeceoul 2999 


Filename: meth_0102.m 
$description: this program uses Methods i and 2 to detect 
various 


$transient signals in a coloree€@ noise environment. 


cme cw mc nr mm wr ww ww eee ee ee ee ww we we we we ie 


clear 

%{s,fs,bps] =wavread(’H:\matiab2\chesis\data\wav_data\record’); 
school 

%{S,fs,bps] =wavread(’C:\My 
Documents\dad\m_files\chesis\data\wav_data\record’); home 

$ioada ‘H:\matlab2\thesis\data\mat_data\fan.mat’ 

load ‘H:\matlab2\thesis\data\mat_data\fan2.mat’ 

Sload ‘C:\My Documents\Gad\m_fiies\thesis\data\wav_data\fan.mat’ 
load ‘C:\Myv Documents\dad\m_files\thesis\data\wav_Cata\fan2.mat’ 
£S=10025; telephone quality 


Sigs eae ze, Ene conridence vector: 
narrowband signais: 

Mio keteormwesees 5 e- oilh=(50;.1:.5]; 
Mogae; 2 |e bzm=(200; .2; .2]; 
nosSovore 05) -mosh=(600; .1:.08]; 
not woe 203 | -mb4h=(900;.1; .05)]-; 
MoS 2.01) nbeosn=(4000. 1-21); 
Moots 000>.1-.01)] ;nb6éh=(5000; -1;.003),; 
NEB=(Mibt nbz -nb3,nb4,nb5,nb6] ; 

New (mela moezhn,nb3h,nb4h, nbosh, nbs6h] ; 
Bee= hee or. 

BBW=[.1;.01]; 

tYials=lengthn (NB); 

Max at r=max (NB’ ); 
Mexeoecr— = Oor (iS max dtr(Z)) ; 


se alef= 
STRSOUND=[] ; 
%CSOUND=[1; 


SWSOUND=[]; 
peewee c=. trials: 
e= 6" 


Grer=Ne (eee): 
£tr=NBH(1,c); 


SaQcranNo 6.6): 
Gtr=NBHC2-e): 


Satr=BSc (li): Sduration Of Colored 26-Se Eranslrent 
Sdtr=BBW(1); Eduiration ©2 woalce nolse transiene 


Saer=Ns(5,c) > 
atr=NBH(3,c); 


Moese Bee (2) 
$atr=BBW (2) ; 


$create transient signal: 

%A) narrowband sinusoidai transient 

NS_T=floor(fs*dtr); t%tnumber of sampies in che transient 
tot_time=NS_T/fs; 

t=0: (cot veime a Nees) ) cote cine. 

V=COS (2"pi* Pere eae 

ham_win=hamming(length(t)); 

y=y. *ham_win’ ; 


%B) broadband transient signal: 

SNS _T=floor(fs*dtr); %*number of samples in the transient 
$yv=randn(1,NS_T); 

$ham_win=hamming (NS_T); 

Sy=y.*ham_win’; 


SC) went ss rome io: 

%NS_T=floor(fs7dtr); tnumber of sampies in the transient 
$corc_time=NS_T/fs; 

he=Ue Eon. time/ (NS 7-1) ) -corEmeime- 

Bub=— (0 Ger, 4) dtry72 3"dtr,/ 2 diri>; & time) breakpoines 

%fp=[0 200 100 150 3003; % instantaneous frequency 
breakpoints 

%p=polyfic(tp, tp,4); % fit 4th order polynomial over time 
Sy COD eee) ; 


%D) colored noise transient signal: 

$NS_T=floor(fs*dtr); tnumber of samples in the transient 
$y=fan2(20000:20000+NS_T-1) ’; 

Sham_win=hamming (NS_T) ; 

S$y=y. *ham_win’; 


%E) colored noise transient signal: ELViS 

tLOade (Ca \My 
Documents\dad\m_files\thesis\data\wav_data\eivis_a.mat’ 
$load ‘C:\My 
Documents\dad\m_fiies\thesis\data\wav_data\elvis_b.mac’ 
load ’H:\matlab2\thesis\data\mat_data\elvis_a.mat’ 
$load ‘’H:\matlab2\thesis\data\mat_data\elvis_b.mat’ 
SNS_T=length(elvis_a); 

$y=elvis_a’; 

$ham_win=hamming (NS_T) ; 

$y=y.*ham_win’; 

talc=o. 


VE=Tane- 
vi (LOOOL (10000; NSet) )=ve(lO00l: (10000+NS _7)) tates: 


%F) use the data fan (includes the tapping) as the input 
Signal 


tyf=fan; 


#mechod 1: pre-whitening 


IZ 


meLZ0nz..l a,zih_a,z21_a,z2n.2,27234_a, z3be oeeeme meee: 01 Filic(y:): 


SleZ Osea, Ziti dg 22a, 220d, 254 vey 73h Ma SNReL=mecrmeel0l fils 
Ay) 3 

SCO ele cin va, Zo ld, 22a. Zoe 7 3n avT=eree Serv.) 
[STG1,STG2,STG3,STG4) =wavepackz (yf) ; 
$SNR_1=my_SNR1 (STG1, STG2, STG3,STG4) ; 

%SNRi=my_SNRi4plots(STG1, STG2,STG3,STG4); 


method 2: post-whitening (Cristi’s algl) 
(Pipes nb, 221. b, 22h. b, zo) _b, zsh =methed 02 filt(yt):; 


peor 2 no, 22) = b,22n. 5,23! Db, Zane bDZSNR Zl =metnod_02222 itv (y 
a be 
%SNR2=my_SNR24plots(y£); 


$method 3: SVD model (Cristi’s alg2) 
ae CN 3 C7 oe Oo CeO C= meenoGs 03. Ever 


éy=(y,zeros(1,max_dtr-iength(y))]; 
$TRSOUND=[TRSOUND;y]; 
SCSOUND=[CSOUND;yf£‘}; 
SWSOUND= [WSOUND; z0‘J; 


SODLO0t resulcs: 

S$3_a=3;13_a=3+8S3_€a; 
e390 D=4-13. b=3+53_)b; 
S4_ a=5;14 a=7+s4_a; 
s4 5=6;14 b=/+s4_bp; 


Evg-plee methol02 (STGl, STG2 (1, - jo STG2(2 2 )] StGs (S320, ) 7 S165 
2) eee es eee 
STG4(s4. a3 ie sle4i{s4 Jb, Set nee 


SE2G=p21 Ot mer ro102 (20) 241. 6,721 nea, 22 ie oy ee ed, Zone, Nowe 
ees pc 
Pig=ploeemetn0l02 (STG, Z21l DrZin ba22) b>, 22h 2b, come, NSS 
Ji pea o ae 

efi G=plot meth03.(z0, 21 eeg2acs7emeg we ewer -Oee, osly eG); 
Progune (216). 

specgram(yf,512,f£s,512,400) ; 

fig=figtl1; 

$end 


4. Plotting Routines 


$Name: John D. Stevens 

Date: 26 Feb 2000 

$Filename: plot_meth0102.m 

SAESCYIDEION: this Fumectr0on will yeigea Ssicma errougmeerte LIVrst 
A2GOrlenam 

that was developed by Dr. Roberto Cristi and pilot it 


% 
% 

D=ploeemecnU 102 (2-0, 22r0c, ZEl parte ee, 2ilo Zeer, NS 1 fo um) 
%5 
ee 
LUunee Lon 


D=pVOESMeCen CMO (chain, 2h. 20.0 Cee 2 heer eee el Se CoCr Undy a 
TTL=( ie 

Ssip-i 

cumilie mranoney= 2 5 


$indexes due to downsampling: 
ind1=10000; 

imndg2=f£ leoortina ly 2): 
ind4=floor(indl/4) ; 
indgdS=floor{indl/s) ; 


Ns_tl1=NS_T; 

Ns@e2—-f£loor(NS-1/2)- 
NS Vt4—fleor (NS 8/4). 
Ns_t8=floor (NS_T/8) ; 


create the spikes surrounding the signal 
Sfar—-Zeros (s1ze(Ztas)); 

Seas ( iG ds)—masxcni2tar) > 

Stan vidal spleen (Ns tl~slip) )=max(ztean): 


sf0=zeros (size (ze£O)): 
S£O(andZ)=max (220): 
SEO (ind2+f locr (Ns.t2* slp) )=maxtZe0) 


s£Ot=zeros(size(zf0t)); 
sft0t (ind2) =max(zi0C) ; 
SfOr(indZz+ELoOor(Nsetz.s! >) jomax (zECE). 


sfl=zeros(size(zfl1)); 
stilt ind?d)=max(z£l):; 
Sftltind4+ilcor(Nsrcea* slp) j=max(ztl)-; 


sflit=zeros (size(zrtle)?): 
sfile{(ind4 jemax(2r bie 
sflt(ind4+floor(Ns_t4*slp))=max(zft1t) ; 


Sfi2=zeros (size(zE2Z) jj; 
sfi2 (ings) =max Uz ic 
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See indo loor (Ns _t8*slp))=max(ze£2) ; 


sf2t=zeros(size(zf2t)); 
Srculineo) —max (| Zi2t); 
SiZetindsS+tloor(Ns_t8*silp) )=max(Zzr2c) ; 


create the time vectors: 
Nar=length(zfar); 
NO=lengEn( Zed) - 
Nl=length(zfl1l); 
N2=length(z£2); 


ES= 2025); 
f_time=Nar/fs; 


r=0:f_time/ (Nar-1):f£_time; 
=0:f time/ (NO-1):£ time; 
=0:f£_time/(N1-1) :£_time; 
=0:f_time/ (N2-1):f£_time; 


plot the figures 

figure (fig_num) 
Pecmucesar,zfar,‘b’,t_lar,sfar, x.) - 

Greeley irr, “AR WHITENING FILTER’ |); 

ylabel (’MAGNITUDE’); 

axist(h0;£ time,min(zfar),max(zftar)]):; 
xlabel(’TIME (seconds) ’); 

S$ylabel (’MAGNITUDE’); 

sLeqenat |’ SNR =." ,numZ2str(SNR(1)),* €aB’))- 


figure (fig_num+1) 

Sopot (274, 1), Dob (tO, 20, bo’ ,t 0, S20, 1-7): 
Stitle([TTL, ‘LOWPASS FILTER STAGE 

1 iere xis O10, f time, min(z£0), max(zre)))— 
title([TTL, ‘LOWPASS FILTER 

(APPROXIMATION) ‘]),axis([0,f_time,0,1]); 

$ylabel (‘MAGNITUDE’) ; 

SLegenadul SNR = “,nUMZSEr(SNR(2Z)),’ dB*i}-: 

Suppl Obie. ,2) ,pLot{t_.0,ZE0t, paeee st OC 2" er 
$title(’HIGHPASS FILTER STAGE 

Wj axis( (0, £.cime,mintzfe0t);max(zeec). ): 

title([TTL, ’HIGHPASS FILTER (DETAIL)’]),axis([0,£f_time,0,1]) 
xlabel(’TIME (seconds) '); 

$ylabel (‘MAGNITUDE’) ; 

Brecgena (|, sNhk = ~,mumZSEV (SNR (3) ). Ges 


figure (fig_num+2) 

SUDpLOEW2 rly lye lOt( Gall, Zr leeeo. Galea i tome 
Stitle([TTL, ‘LOWPASS FILTER STAGE 

2 pepe ccroug i EW CaAme. mae ey) Max 2 Ll) yi 
title([TTL, ‘LOWPASS FILTER 
(APPROMIMAT ION) ))> axis (.[0,2 2erme.e... |: 

%ylabel (‘MAGNITUDE’) ; 

legend ((’ShR-=."j,numZ2str(SNRi4)),° Ge’ )); 

Su plCrt ieee me DlOE( GC. 1) 2il te. < ee Lior tes ee 
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$title(’HIGHPASS FILTER STAGE 

2, Axis (00 seeine, Minter lc) mesa ae 
title([TTL, ’HIGHPASS FILTER (DETAIL) 
xlabel(’TIME (seconds) ’); 

$ylabel (‘MAGNITUDE’); 

Slegend( |’ SNR =" *,nume@ser(SNR(5)), ‘oda Ie 


ree al 
a een ta! tf) tame, 0,1 )4-> 


figure(fig_num+3) 

SUbDIOEZ ~!) Sloe Cem 2i2, DY eae ems roa 
S$title((TTL, ’LOWPASS FILTER STAGE 

3 | )peNeSs (10. fe ce imiewme a 262) anata oe 
title([TTL, ‘LOWPASS FILTER 

(APPROX® MATION) ) \irexist (Ot time, Onli 

%ylabel (‘MAGNITUDE’); 

Sslegend([‘SNR = “,num2str(SNR(6)),° aB’]); 
subplot(2, 1,295, DloOu(t. 2, Z2iZtiy + pte sh eee: 
$title(’HIGHPASS FILTER STAGE 

Ss )waxis (04 ft came min(Ze2c) max (ees e 
title( (TTL, ’HIGHPASS FILTER (DETAIL) ’ 
xXxlabel(’TIME (seconds) ‘); 

tylabel (‘MAGNITUDE’); 

$iegend([’SNR = ‘’,num2str(SNR(7)),‘’ aB‘’]); 


a 
Vy, acre GV hOtet ec ame, 0, 1.) ) 


b=fig mums: 


¢Name: John D. Stevens 

mDates 22 cue 999 

S$Filename: SNR_plots3.m 

description: plots the results from my_SNR1 and my_SNR2 


clear 

S$load up results: 

load ‘H:\matlab2\thesis\data\mat_data\SNR_ARb.mat’ 
load ‘H:\matlab2\thesis\data\mat_data\SNR_Mlb.mat’ 
load ‘H:\matlab2\thesis\data\mat_data\SNR_M2b.mat’ 

Bu Oe ee: \My 
Documents\dad\m_files\thesis\data\wav_data\snr_ar.mat’ 
load ‘C:\My 
Documents\dad\m_files\thesis\data\wav_data\snr_Ml.mat’ 
load ’C:\My 
Documents\dad\m_files\thesis\data\wav_data\snr_M2.mact’ 


FITR-O; 102 3000. 

ATR=0222002 1; 

thresh=3; threshold in dB 
ax=[min(FTR),max(FTR),min(ATR) ,max(ATR),-20,50]; 


figure (1) 

surf (FTR,ATR, SNR_AR) ;axis(ax); 
mesh(FTR,ATR, SNR_AR, SNR_M1);axis (ax); 

C=colormap; 

title(’NARROWBAND DETECTION REGION FOR WHITENING FILTER’); 
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xlabel(’FREQUENCY OF TRANSIENT (Hz)'); 
ylabel (’AMPLITUDE OF TRANSIENT’); 
Zlabel(’SNR (dB) ‘); 


Eilgure (2) 

$SurI(FTR,ATR, SNR_M1);:axis(ax); 

mesh(FTR,ATR, SNR_M1,SNR_M1) ;axis(ax); 
title(’NARROWBAND DETECTION USING WAVELET DENOISING’):; 
xlabel(’FREQUENCY OF TRANSIENT (Hz)’); 
ylabel(’'AMPLITUDE OF TRANSIENT’); 

zlabel(’SNR (dB) '‘); 


figure (3) 

Ssurf(FTR,ATR, SNR_M2) ;axis(ax); 

mesh(FTR,ATR, SNR_M2,SNR_M1) ;axis(ax) ; 

title(’NARROWBAND DETECTION USING MULTIRESOLUTION INNOVATION - 
THRESHOLD = 3 dB’); 

xlabel(’FREQUENCY OF TRANSIENT (Ez)’); 

yYlabel( ‘AMPLITUDE OF TRANSIENT’); 

Zlabel(’SNR (dB) ’); 


figure (4) 

X=SNR_AR; 

iE idiciccc thresh) 
x(I)=thresh; 
Gomtourr (FIR, ATR-x);colorbar: 
GOlonmap (jet); 
BOeGOLOr (PT R,ATR x): 
title(’NARROWBAND DETECTION REGION FOR WHITENING FILTER - 
THRESHOLD = 3 dB‘); 
xlabel(’FREQUENCY (Hz)’'); 
ylabel (‘AMPLITUDE’); 


figure(5) 

Solin eM: 

iL—pes jaye (be jelebgercisy) | 
x(1)=thresh; 

GOmMecourr (PTRAATR, x) colorbar;: 
colormap (jet) ; 
meGOlOr (PIR, ATR, x); 
title(’NARROWBAND DETECTION USING WAVELET DENOISING - THRESHOLD = 
Swag): 

xlabel(’FREQUENCY (Hz) ’‘); 
ylabel (‘AMPLITUDE’); 


figure (6) 

SIS OS 

P=find(x<thresh)- 

x(I)=thresh; 

Concours: (FTIR, Als, x); cOlLOmoaT 

colormap (jet); 

£0COlOr(PIRAATE =) 

title(’NARROWBAND DETECTION USING MULTIRESOLUTION INNOVATION - 
THRESHOLD = 3 <8’); 

xlabel (‘FREQUENCY (Hz)’); 


ey 


ylabel (‘AMPLITUDE’); 


SDETERMINE PERCENT IMPROVEMENT ON DETECTION: 

Sea) FOR THE ENUERE SoBCrTRUM: 

tot_pts_AR=prod(size(SNR_AR) ); 
det_pts_AR=find(SNR_AR>=thresh); find all points >= 3db; 
det_pts_AR=length(det_pts_AR) ; 
pct_det_AR=(det_pts_AR/tot_pts_AR)*100; %tpercent detection for 
whitening filter 


tot_pts_Ml=prod(size(SNR_M1) ); 

det pts =Ml=find (SNReM@—cehresh): ¢findvall” points >=) 3do. 
det_pts_Ml=length(det_pts_M1); 

pcE_ det_Ml=(det prs MI /tot pts Mi) “0G es -sercene detecevony cor 
whitening filter 

pet_incrM1=100* (pct_det_Ml-pct_det_AR)/pct_det_AR; %tdetection 
percent increase 


tot_pts_M2=prod(size(SNR_M2)); 
det_pts_M2=find(SNR_M2>=thresh); %find ali points >= 3db; 
det_pts_M2=length(det_pts_M2); 
pcet_det_M2=(det_pts_M2/tot_pts_M2)*100;: %percent detection for 
whitening filter 

pct_incrM2=100* (pcet_det_M2-pct_det_AR)/pct_det_AR; %detection 
percent increase 


POTlODET “TOT [| pecracec An; pet det _Ml-= pet det M2] 
Pe? Bee R= (0; peru neu pet incrM2 | 


$B) FOR LOW FREQUENCIES BELOW 2kHZ ONLY: 
fin=find(FTR==2000); %tcolumn where freq = 2000 Hz 


POEUZ Kk AR=preod(stzZeUSNR AR. be fini) 

@eta2k AR=find(GNR AR(-, 1: fin)S=thresh);> —~tiind all 'peines >= 
Selon 

det_2k_AR=length(det_2k_AR) ; 
pct_2k_AR=(det_2k_AR/tot_2k_AR)*100; %*percent detection for 
whitening filter 


EGuscheMl—prod(sizet Shinai (:, l:fin) 

deteck Mi=ftind(SNRIMOG | fie thresh) ctind all coints >= 
Sisley 

det_2k_Ml=length(det_2k_M1) ; 
DCE_2ZK_Ml=(det 2k MIl/tet_2k Ml)*100;> tpoercene cetection for 
whitening filter 

pet_2k_incrM1=100* (pet_2k_Ml-pct_2k_AR)/pct_2k_AR; %detection 
percent increase 


EOE 22M? =pred (size (Shine rr) ) 

det.2k M2=find(SNR_M2 (32 )-fin)s=thresn) = $find ait  eointes == 
erelor, 

det Z2ZkaM2=length (det 22a MZ) 

pct_2k M2=(deEel2k M2/tee 2k M2) *100; tpercene detection for 
whitening filter 
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PGbez kKeimerM2=100* (pcet_2k_ M2-pct.2k_AR) /pct_2k AR] sdecection 
percent increase 


PetepEtT Zk=(pect_2k AR-pct_2k_M1;pct_2k M2} 
PeteSe wile RezK=|O7pct_2k.inerMl- pct 2k ineri2 


$C) FOR LOW FREQUENCIES BELOW 300 Hz ONLY: 
fin=find(FTR==300); %column where freq = 300 Hz 


tot_3c_AR=prod(size(SNR_AR(:,1:fin))); 
Germoecwer-ringd | SNR. AR(:;,12¢fin)s=thresh);  ‘sfind ali poinss >= 
BSS). 

det_3c_AR=length(det_3c_AR) ; 
pet_3c_AR=(det_3c_AR/tot_3c_AR)*100; %percent detection for 
whitening filter 


POtresem))l=prod(sizetSNR_M1 ( ee) 

feeeac MIl=tinag (SNR _Ml( =) 1 fin) ==thresh) er tion, 2! points >= 
Seb- 

det se Mi=length(det_3¢_ M1); 

Petese Mi=(cdet sc M1/tot_3c_Ml)*1l00- tpercent adeteceron seer 
whitening filter 

Bet seprmerMI=100* (pct serMl=pct seo Ab 7 pet Seen, *s0Ceceezi1o) 
percent increase 


tOumssGuM2=prodad(size(SNR_M2(: 1. fan) ) je 

det serM2=find(SNR_M2(:,1:fin)s=thresh)s ~sei0d a.) poles >= 
LXelley- 

eepese M2-length(det_3¢_M2) ; 

Peeece M2Z2—-(det_3c_M2/tot_3c_M2)*100; tpercent Gdecectioy for 
whitening filter 

WewmecmeimerM2=100* (pct_3¢_M2-pets3c_AR)/ pe a Cuan ce Co.) 
percent increase 


PQIm@UE tase= (OC ese AR; pCt-scC_ Mil; pct s3e uz) 
PCT BETTER. 3c=[0;pct_3¢_ incu peewee ae | 


Se 


Di Display Filter 


$Name: Jorn D. Stevens 
Date: 2i Jan 2000 


Filename: smoother.m 
$description: this program will take a white noise process and 
Gusout 


$a normalized smoothed version for display. This process aiso 
$gives an approximation to the noise variance and the signal 
variance, 

$which is useful to compute the signai to noise rati 


% zf=smoother (y) 


FUNGE LON 62 ir -SmOoOEner (ye 
if Size(y,2) >size(y;1). Y=) > end % make y into column vectors 


Ze ..072; 
=a 
A=[1,-.99]; 


ZE=Piel Cem wel 99 lets 

zf=zf'; make output in row vectors 
$normalize the resuits: 

MZ =Meec( 2 © alee)? 

Zh=7F e/MZE (ones (lengeimtezry 1): 

St 2ey tr obi | tee er Behe) ees. 


6. Broadband Transient SNR Test 


‘Name: John D. Stevens 
$Date: 22 Jul 1999 
$Filename: SNR_elvis 
$description: 


ccm cm cr cr ccm crm we ee es 


clear 

load ‘H:\matlab2\thesis\data\mat_data\fan2.mat’ 

load ‘H:\matlab2\thesis\data\mat_data\elvis_a.mat’ 

load ’C:\My Documents\dad\m_files\thesis\data\wav_data\ian2.mat’ 
tload ’C:\My 
Documents\dad\m_files\thesis\data\wav_data\elvis_a.mat’ 


fs=11025; S$teLiephone quality 
SNR_AR=[]; 

SNR_M1L=[]; 

SiR z=); 

acr=)1 
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$NS_T=length(elvis_ a4); 
NS_T=floor(fis*dtr); %number of samples in the transient 
ham_win=hamming(NS_T) ; 


SELVIS 

$y=elvis_a'; 

Sy=y/max (y) ; Snormalize the elvis signal 
$y=yv.*ham_win’; 


SWHITE NOISE: 
y=randn(1,NS_T); 
y=y.*ham_win’; 


fom atr=0:.005:1; tamplitude 
eve 
add transient to benchmerk signal 
Vie= mane 
me(LOOOL:{ LGM00+NS T))=yiel0001: (LOQO00NS_T) )eatr*y- 


method 1: pre-whitening 
(STG1, STG2, STG3 , STG4] =wavepack (yf); 
snr_l=my_SNR1(STG1,STG2,STG3,STG4) ; 


method 2: post-whitening 

Bie e2=my eSNR2Z (ye) - 

SliReeR=(SNR_AR:>Snr_1(1)]; 
SNReOML=(SNR_-M1;snr_i1(2) ]; 
SNR_M2=[(SNR_M2;snr_2]; 

end 

Pea Orria 5.) * 4 

res=[SNR_AR, SNR_M1,SNR_M2]; 

figure (16) 

plot (ATR,res,ATR,3*ones(size(ATR)),‘k:'); 

axis([0 1 -20 max(SNR_M1)]); 

Siebel iP LITUDE OF ELVIS”); 

ylabel(’SNR (dB)‘); 

Megend(['’Whicening Filter’) 77 Method 1’), [’Method Zayy. 
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B. 


METHOD 1 


1. Method 1 Filter Process 


Name: John D. Stevens 
$Date: i March 2000 


$Filename: wavevack.m 
$description: this program will first whiten a colored noise 
Signal 


y and then decompose it into a Haar wavelet decomposition packet 
$framework the channels will be white noise processes 


function [STAGE1, STAGE2, STAGE3 , STAGE4] =wavepack(y) 

if size(y,2)>size(y,1), y=y’;end%make the signal a column vector 
perform ARMA process to get a whitened signal 
[b,a]=cristi_ar(y); 

eQ=filter(a,b,y)’; 


@filter the signal through the Haar wavelet transform packet: 
Mol i/saqrt (2) sare (2) ly 1 /sqree meee SOrc (7 ae 
temp=M*reshape(e0,2,length(e0) /2); 
Al=temp(1,:); 

Dl=temp(2,:); 
temp=M*reshape(Al,2,length(Al)/2); 
RAD =Cempi l=) 

DA2=temp(2,:); 

temp=M* reshape (AA2,2,length(AA2) /2); 
AAA3Z3=temp(1,:); 

DAA3=temp(2,:); 

temp=M* reshape (DA2,2,length(DA2)/2); 
ADA3=temp(1,:); 

DDA3S=temp (2, :); 

Femp=M" resnape (bl, 2, tengChi(Dl)7 2); 
ADZ=temp (1.2); 

DD2=temp(2,:); 

temp=M* reshape (AD2,2,length(AD2)/2); 
AAD3=temp(1,:); 

DAbpS—cemp (2, :)- 
temp=M*reshape(DD2,2,length(DD2)/2) ; 
ADD3=temp(1,:); 

DDD3=temp(2,:); 


STAGE1=e0 ; 

STAGE2=[A1;D1]; 

STAGE3=[AA2;DA2;AD2;DD2]; 

STAGE4= [AAA3 ; DAA3 ; ADA3 ; DDA3 ; AAD3 ; DAD3; ADD3;DDD3] ; 
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$Name: John D. Stevens 

$Date: i March 2000 

$Filename: wavepackz.m 

description: this program will first whiten a colored noise 
Signal 

$y and then decompose it into a Haar wavele= decomposition pacxet 
$framework. The result will include the display filtering 
process, 

$smoocther 


function [STAGE1, STAGE2,STAGE3, STAGE4] =wavepack (y) 

if size(y,2)>size(y,1), y=y’;end%émake the signal a columm vector 
toerform ARMA process to get a whitened signal 

Pome cristi arly); 

e@=fialter(a,b,y)’: 


filter the signal through the Haar wavelet transform packet: 
Medeecemet 2), L/saqrt(2)-1/l1/saqrt (2), -1/squcl(2- 
temp=M*reshape(e0,2,length(e0) /2); 
A femp ge. 2) 

Di-cemoa(2, :); 
temp=M*reshape(Al,2,length(Al)/2); 
AA2=temp(1,:); 

DAZ=temp(2,:); 

temp=M* reshape (AA2,2, length (AA2) /2) ; 
AAA3=temp(1,:); 

DAA3S=temp(2,:); 
temp=M*reshape(DA2,2,length(DA2)/2); 
ADA3=temp(1,:); 

DBAS =temp (2,2 )'; 
temp=M*reshape(D1,2,length(D1)/2) ; 
Ze =temoil,:)' > 

DD2Z2=temp (2; :); 

temp=M* reshape (AD2,2,length(AD2) /2) ; 
ARDS =eemp (lye) : 

DAD3=temp(2,:); 

temp=M* reshape (DD2,2,length(DD2) /2); 
ADDS =cemp (i); 

DDD3=temp(2,:); 


STAGE1=smoother(eQ) ; 

STAGE2=smoother([A1;D1]); 

SiGe s=smoocmer ( (AA2;DA2;AD2;DD2]); 

STAGE4=smoother ( [AAA3 ; DAA3 ; ADA3 ; DDA3 ; AAD3 ; DAD3; ADD3; DDD3] ) ; 


2: Narrowband Transient SNR Test 


$*Name: John D. Stevens 
SDacte? fl Jan 2ZO0Cc0 
Filename: my_SNR1.m 


Saescription: find the s2qnal Eo noise ratio of the all ere 
channels 

$in the Haar Wavelet Packet 

% SNR=my_SNR1 (STAGE, STAGE2, STAGE3 , STAGE4) ; 

% 

NN Ned oN me sf me ln Gael fp nde |) Se ea rnc eed ee fe fee eel ede 


function SNR=my_SNR1 (STAGE1, STAGE2, STAGE3 , STAGE4) ; 

Sthe assumption is that the transient signal is 

NS. T=floor (2 i an: 

slp=1; 

NSite lisa. 

Ns_e2=t loom(NsS 272). 

Ns_t4=floor(NS_T/4) ; 

NS JEG=f loom Noe, oc) ; 

$zfar 1s a column vector while all the other inputs are row 

vectors 

$indexes due to downsampling: 

ined! —100 Co, 

Tnddiwrin—Ind lst loer (NS etl *s tae 

IndZ foo, (indi 2); 

indZeeimei nds +E loor (NS e2*sipr 

Imada moor (imal 4): 

inde) fim—-ind4+tloor (Nsee4-sip) 

indé=floor (ind, 3); 

INGS  fin=itnast+Eloor (Ns cSs*sia 

S=STAGE1; 

Al=STAGEZ (1,:); 

DIi=STAGEZ (2). 

RA2=STAGE3 (1, 

DA2=STAGE3 (2, 

AD2=STAGE3 (3, 

DD2=STAGE3 (4, 
( 
( 
( 
( 
( 
( 
( 
( 


as ™s 


seconds long 


Ses 
eS, 
one 
ee 
AAA3=STAGE4 
DAA3=STAGE4 
ADA3=STAGE4 
DDA3=STAGE4 
AAD3=STAGE4 
DAD3=STAGE4 
ADD3=STAGE4 
DDD3 =STAGE4 
k= 1S; 

nS=var ([ 
nAl=var ( 
nDl=var ( 
nAA2=var 
nDA2=var 
nAD2=var 
nDD2=var 


se 


™=o 


. 
Fd . 


« moe 


“=e 


e 
a s 


=e 


e 
F . 


=e 


e 
Fd e 


oe 
Seer ee” ae See eel bal al Saal 
= 


jy 
Ze 
5 
ee 
5 
6 
a 
8 


™=e 


s . 


ire moan Loom ia" ind lesran).- bengtint Ss) pili) + 

11nd Se et tloor{ k*tnmd2 fin): length Al): ]): 
1 ( Si aaa) ae Mengteh (D1) ) 1); 
2(1: (ind4-1 
Oe ab ave ks! al 
(1: (ind4-1 
Calee st 


] 


2(tloor(k*ind4 fin) = lengeniAAZ)) 
DA2 (Eloor (k+ind4 fin) Length (DA2) 
WADZ(E1lOOr( 7 1ne4 erin) veng en (Ab2)) 
an pee i rar cree eenae 


DA 
AD ; 
DD 


) 
ee 
) 
) 


ind4-1 


e 
s 


ot 
[A 
[D 
(TAA 
([ 
([ 
al 


)), 
Ee 
)) 
)) 
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nAAA3=var ([AAA3 (1: (ind8- 

i pa nomr LOOT (K*inds Ein): bengen(AAA3))]); 
nDAA3=var ({[DAA3 (1: (ind8- 

ip, DAA SX rloor (k*inds fin): length (DAA3) )])}:; 
nADA3=var ([ADA3 (1: (ind8- 

ee SPAS tloor(k*ind8’ fin) :lengtn(ADA3))]); 
nBPAS=var( {DDA3S (1: (inds-— 

Pees OOor (k*inds fin) -léngeh(DDA3)) jo 
nNAAD3=var ({AAD3 (1: (ind8- 

ip SADS( tloor(k*inds fin) :length(AAD3) ))); 
nDAD3=var ({DAD3 (1: (ind8- 

ioewAbSs (tf loortk*inds fim = engten(DAD3) ))]); 
nADD3=var ({[ADD3 (1: (ind8- 

itPrmaADbs (tT loor (k*ands fim) length (Abb3) ))] ); 
nDDD3=var ({[DDD3 (1: (ind8s- . 

preemies (floor 1906 fin) +Lengten(DDD3) ) ]); 
%FIND THE SIGNAL POWER: 
Ss-Vom(o(indl ind] fimo 
sAl=var (Al (ind2:ind2_ fin)); 
Shan (Dl (imdZ2 ame. fin) )- 
SAA2=var (AA2 (ind4:1nd4_fin)); 

SDA2=var(DA2 (ind4:ind4_fin)); 

SAD2Z=vVar (AD2(ind4:ind4 fin)): 

SpMoz—var(Dp2 (ind4-ind4 fin) js 
sSAAA3=var (AAA3 (ind8:ind8_fin) ); 
SDAA3=var (DAA3 (1nd8:ind8_fin)); 
SADA3=var (ADA3 (ind8:ind8_fin)); 
SUDAs=Var(DDAS (1nds8:ind8._fin)); 
sAAD3=var (AAD3 (ind8:ind8_fin)); 

SPAS —=Var (DADS (1nd8:inds fin)}- 
SADD3S=var (ADD3 (ind8:ind8_fin)); 
SDDD3=vVar (DDD3 (ind’:ind8_fin))-; 

get the raw ratio: 

snr_S=sS/nsS; 

sow Al=sAl/nAl; 

SrnieeDil=<c D7 nis: 

snr_AA2=sAA2/nAA2; 

snr_DA2=sDA2/nDA2; 

snr_AD2=sAD2/nAD2; 

Smee =sDU4/mDDZ- 

snr_AAA3=sAAA3 /nAAA3; 

snr_DAA3=sDAA3/nDAA3; 

snr_ADA3=sADA3/nADA3; 

snr_DDA3=sSDDA3 /nDDA3; 

Sat efAD3=sAADS /nAAD3; 

snr_DAD3=SDAD3/nDAD3; 

snr_ADD3=sADD3 /nADD3; 

snr_DDD3=sDDD3 /nDDD3; 
c=[snr_S;snr_Al;snr_D1;snr_AA2;snr_DA2;snr_AD2;snr_DD2;... 
snr_AAA3;snr_DAA3; snr_ADA3;snr_DDA3; snr_AAD3; snr_DAD3 ;snr_ADD3;sn 
Pawo DS |%; 

X=C;X=x-ones(Size(xX)); 

f=fpind (x01) > 

so ( 1) = 01l-SNR=10*1logl0 (x); 

SNR=[SNR(1) ;max(SNR(2:15))]; 


METHOD 2 


1. Method 2 Filter Process 


$Name: Jokn D. Stevens 

$Date: 12 Jan 2000 

$6Fiiename: method_02 filt .m. 

SdeScrivtion:) tChissfunce2onywi)) Pusan sitcona: Eerougn EPes ri rset 
Gude = ve im 

S$that was developed by Dr. Roberto Cristi 

% 

(Zi ze, Zz2) 226723 ,23 | =metnecml cues te (y) 


oP oP 


mc mr re es ae ee ia ee i ee 
a 


EUMGCE On wizie- Zz le 22 22075, 250 | =mMernoceO rie 


hE SPZE(V 1) =size(y,2), y=y’ > end % y row vector 
k=1/sqrt (2); 
tk=. oe 


lena t=Crrsti arte) 

[bly al ble ,alt)]=modeliZki(b,ak) 
(D27a2 OZ aZzt | —=modelZk (ol ,aleek) 
[634435 5082,a3t | =modelZkib2,az, ki: 


Mae le el 
yl=M*reshape(y,2,length(y)/2); 
V2Z=M*reshapety!l (1, -:),2, Lengen(yl) 72) 
y3=M* reshape (y2(1,:),2,length(y2)/2); 


seQ=filter(a,b,y); 
St=tritcer (al bl yitl ))- 
el@=krleerlalt,blt,yl(22.) )74srow vector 
e7=fiirer(az,b2,y2(1,2))e 

eZ t=pmbbern(azt,b2t,y2(25.)) > 
es=filter(as,b3,y3 (1, :)); 
Ssbariheervast, Ost, ysi27 0: 


$zQ0=smoother (e0); 
zl=smoother (el); 

Z V=2 ey Mase Cz |): 
zlt=smoother(elt) ; 
ZIC=2 ie mes ¢zZ Ve): 
z2=smoother(e2); 
Z2=20 (ase 2a. 
z2t=smoother(e2t) ; 
Z2C=ZZt/max(2Ze)- 
z3=smoother(e3); 
Z23=23/max(z3); 
z3t=smoother(e3t); 
z3t=z3Umax(z3t); 
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2. Multiresolution Innovation Filter Bank Routine 


‘Name: John D. Stevens 

$Date: 23 Feb 2000 

$Filename: model2k.m 

$description: This program appiies the multiresolution 
innovation process 


EUNCGEVOMUMmLO2,d2,o02t,azt | =moderZkino,a,k) 
% [b2,a2,b2t,a2t] =model2(b,a) 
* "t" Stands for "télde”, whrenh 1s the high pass, 22lce- 


[Ce, Genk zil=peivd(b) a); 
if length(Ce)<length(Co), num=k*Co+k* [Ce, zeros(1,length(Co) - 
length(Ce))]; end 
1f length(Ce)>length(Co), num=k*Ce+k* (Co, zeros(i, length (Ce) - 
length(Co))]; end 
1£ length(Ce)==length(Co), num=k*Ce+k*Co; end 

den=Fz; 
if length(num)>length(Fz), den=[Fz,zeros(1,length(num) - 
length(Fz))]; end 

[F1,G1,C1,D1])=tf2ss(num, den); 


if length(Co)+l<length(Ce), num=k*Ce+k*[(0,Co,zeros(1,length (Ce) - 
length(Co)-1)]; end 
if length(Co)+1>length (Ce), 


num=k* [0,Co])+k*[Ce, zeros(1,length(Co)+l-length(Ce))]; end 
if length(Co)+1l==length(Ce), num=k*Ce+k*[0,Co]; end 
den=Fz; 


if length(num)>length(Fz), den=[den, zeros(1, length (num) - 
length(den))]; end 
if length(num)<length(Fz), num=[num, zeros(1,length(den) - 
length(num))]; end 

[P2,6G2,C2,D2]=tti2ss (num, den); 


% overali state space 

Meelengem(—2)=-lengen(Pl) +1, 
Fi=([Fl;zeros(1,length(F1)-1),1],zeros(length(F1)+1,1)]; 
Gi= ClO]; 
Ci (Cis 0]: 

end 


A=F1'; 

Baeveé2 ,Cl i]; 
C=G2’'; 

D= (Don) :: 
Q=b"B- 
R=D*D - 
S=B*D’; 


[K,P] = mydige(A, eye(size(A)),C,Q,R,S); 
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Pez, a2l=ss2ti (45k 7c ee 
Sagma—-sqre (C7 Fe Fae 
b2=b2/sigma; 


$ "tilde" - the high pass filtering 
af length(Ce)<length(Ce),) num=k*Ce=k* (Ce7zeros( i, lengen (ee) = 
length(Ce))}]; end 
if length(Ce)>length(Co), num=k*Ce-k* [Co, zeros (1, length(Ce) - 
length(Co))}]; end 
if length(Ce)==length(Co), num=k*Ce-k*Co; end 

den=FzZ; 
if length(num)>length(Fz), den=[Fz,zeros(1,length(num) - 
length(Fz) ))> end 

[Pip Gl Ci Di \=eezee aum, den) 


if length(Co)+l<length(Ce), num=k*Ce-k*[0,Co, zeros(1,length(Ce)- 
length(Co)-1)]; end 
1f length(Co)+1>length(Ce), num=k*[0,Co]- 
k* [Ce, zeros(1,length(Co)+1l-length(Ce))]; end 
1f length(Co)+l==length(Ce), num=k*Ce-k*[0,Co]; end 
den=F2Z; 
if length(num)>length(Fz), den=[den,zeros(1l, length (num) - 
length(den))]; end 
if length(num)<length(Fz), num=[num,zeros(1,length(den) - 
length(num))]; end 
([F2,G2,C2,D2]=tf2ss (num, den) ; 


oP 


overall state space 

1f length(F2)==length(F1)+1, 
Fi=((F1;zeros(l,length(Fl1)-1),1],zeros(length(F1l)+1,1)]; 
G1=[G1;0]; 
Ci= hel, Ot; 

end 


A=F1*: 
B=lG2Z. 2 Cl | 
G=CZ. 

D=|sD2 oi 


O=B*B" ; 

ReD~pD* ; 

S=o 2) a, 

[K,P] = mydlqge(A,eve(size(A)),C,0,R,S); 
\b2@aeZztlt=ssl2ert A K,C rl) 


Sigma=sqrt (C*P*C’+R) ; 
b2t=b2t/sigma; 
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3. Narrowband Transient SNR Test 


Name: John D. Stevens 

$Date: 13 Jan 2000 

$Filename: my_SNR2.m 

description: this function wiil determine the SNR of various 
transients in white noise 


ve 


SNR=my_SNR2 (y) 


ae cP 


: ee ee ee ee ee ee et 
function SNR=my_SNR2 (y) 

if size(y,1)>size(y,2), y=y’; end *& y row vector 

kane SGraee |. 

tik= 5 


[Sy alecriStivarty’): 

(bivvad blt alt |=meodelek(byaskk). 
[b2,a2,b02t, a2t])=model2k(bl,al,kk); 
[pb37-e3,b3t, azt l=mocel2k(b2,a2, kK); 


Melek =i ee 1 
yl=M*reshape(y,2,length(y) /2); 
V2="-reshnapet(ye(l,:),2,Llengen(vii7 2 
y3=M* reshape (y2(1,:),2,length(y2)/2); 


e0=filter(a,b,y); 

Cl—tadtter (all, yitl) <<) ); 
elt=fidtterltalt, DIE, vlL(274) )4) srow veccor 
Cvaiibeentas, D2. Veil. s ae 
eve=tilterlatteavit. yel2, .))> 
es=ftalter (as), bs, yo(l, 2) > 
Sst—flliemiastynst, youl. 7) 7s 


$find SNR: 
NS2T=£ loom (ke Oe Sie 
slp=1; 

Nemec l=NSel- 

NS vez—e1oor(Ns T/2); 
Memeo —t Gon (NS21/4); 
Nswec=£loor(NSiT/S) : 


$zfar is a column vector while ali the other inputs are row 
vectors 

Sindexes due to downsampling: 

indl=2#0000; 


inal: fin=andl+?loor (Ns 2b l*sip); 
tmd2 = LOor line ly 245; 
ma. Lin=ind2+t Loor (NS2t2*sip)- 
PHed=f bLOoor (ama. / 4)" 
gnedefin-iniea- floor (NS ta=sip) 
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PWG —f1 OS (ame 6). 
Mico fin=i1nde+rt loor (Ns te ~sip)-? 


ZEQ=eL- 
z£0t=elt; 
zf1l=e2; 
zflt=e2t; 
Zee =63- 
Becta e J 


a 3 Ss 

find the variance of the noise 

Ingvavam (ZEOC(l tind2—li ZE0 (iE leoor (ht reget nn, - lenge (zi). |). 
mivec=var( (zftOt( l=: (2nd2Z— 

Ly), ZEQG( Eoor (kk ind2eurin) length (zr0e een. 

moil=var((zil(l: (ind4-1) ), ZE1( floor (k*indtetin) -tengen(zel)) |); 
mnit=var([(z£1lt(1: (ind4- 

ij cele te loc We ater in) = length (Zee) oe 

MiZ=Var((zi2tl a anee—l)) 7 Z2i2it loon kh inde wtin = lengen (zizZ)))) ; 
MaZtaovar((ZEZze( (inde — 

1) eZee (te loor (k* inde fin) lengen(z1Z2 te ene 


$find the variance of the spike in the region: 
Mo@=vearizZe0 (anagsand2 fin) )> 
mp0G=var (zElt (gnd2-i1nd2 fin) )< 
mpl=var(zf£l(ind4:ind4_fin) ); 
Mplt=svan(zile(indd -indésein) )e 
mpZ=Var(zZrZ(indecinde fan) ) - 
Mp2t=Vainl ZEZe (inae ands fin) ) ; 


eget the raw ratio: 
snr _O0=mp0/7mn0; 
snr_Qt=mpo0Ot/mn0t; 
snr_l=mpi/mni1; 

Sbe dhe—ijete yqedlics. 
snr_2=mp2/mn2Z; 

Sie Zt=moZt/mn2Zt : 


Sonn Snr ar,snr 0 snr Oc ,snr_ |, snr. le vshipes Sue Ze) 
SNR=(Ssnir_0, snr Ct, snr il -snr lt,snr 2, snnezci- 


beemeck £Or-i2il conditions 
x=SNR; 

x=x-ones(size(x)); 

L=Ppined (cee |.) 

oo) =e 

SNR=10 Vogtle (sc)> 


SNR=max (SNR); 
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